diff --git a/README.md b/README.md index 0568295..bfa44e6 100644 --- a/README.md +++ b/README.md @@ -24,9 +24,6 @@ What currently works is - Block Term Decomposition (BTD) - Join Decompositions -See [PIPELINE.md](docs/src/PIPELINE.md) for the current execution flow. - ---- The next updates will include @@ -39,9 +36,21 @@ The next updates will include - Partially Symmetric CP - Tensor Trains -
+See [PIPELINE.md](docs/src/PIPELINE.md) for the current execution flow. + +--- + +## Installation + +Install the current development version through the Julia package manager: + +```julia +add TensorKitchen +``` + +--- ## Canonical Polyadic Decomposition (CPD) Here is how to approximate a tensor `A` by a CPD of rank `r`. diff --git a/docs/PIPELINE/index.html b/docs/PIPELINE/index.html index ab26047..a214121 100644 --- a/docs/PIPELINE/index.html +++ b/docs/PIPELINE/index.html @@ -1,2 +1,2 @@ -Pipeline · Documentation
GitHub

TensorKitchen Pipeline

This document explains how public APIs route into models, solvers, and result converters.

Public entry points

  • cpd(A, r; ...) -> CPDResult
  • nncpd(A, r; ...) -> CPDResult
  • btd(A, blocks, ranks; ...) -> BTDResult
  • tucker(A, ranks; method=...) -> TuckerResult
  • approx(...) -> ApproxResult or auto-routed CPDResult/BTDResult

Default behavior (quick reference)

  • cpd(A, r):
    • init = :alswarm
    • solver = :rgd
  • nncpd(A, r):
    • init = :alswarm
    • solver = :rgd
  • btd(A, blocks, ranks):
    • init = :alswarm
    • warm_steps = 200
    • warm_init = BTDHOSVDMultistartInit(candidates=64, screening_steps=10, block_maxiter=12)
    • warm_rel_error_gate = nothing (run manifold refinement by default; set e.g. 5e-2 to short-circuit on poor warm starts)
    • solver = :rgd
    • final BTD-ALS polish enabled by default for non-ALS solvers
    • max_stagnation_restarts = 1 (retry with stronger multistart when ALS fit-change stalls at high rel-error)
  • tucker(A, ranks):
    • method = :sthosvd
  • approx(model::JoinModel):
    • init = :alswarm
    • warm_steps = 500
    • solver = :rgd

Core execution architecture

Most optimization APIs share this core pattern:

  1. Build a model (JoinModel + backend)
  2. Call _solve_model(...)
  3. Convert to a public result struct

_solve_model lives in src/solvers/solve_dispatch.jl and is the common symbol-to-solver dispatch layer (:rgd, :rcg, :lbfgs, :als, :btd_tsd).

API flows

CPD (cpd, nncpd)

cpd(A, r; ...):

  1. Build JoinModel(A, r; geometry=...) with CPDBackend
  2. Normalize/validate options (solver, geometry, gradient_mode, normalization policy)
  3. Solve through _solve_model(...)
  4. Optionally run nonnegative ALS polishing (for selected nonnegative paths)
  5. Convert to CPDResult

Notes:

  • :als means CP-ALS.
  • Manifold solvers (:rgd, :rgd_fixed, :rcg, :lbfgs) share dispatch with other pipelines.

BTD (btd)

btd(A, blocks, ranks; ...):

  1. Build a uniform Tucker family via TuckerJoin(...)
  2. Wrap as JoinModel with BTDBackend
  3. Choose effective initializer:
    • solver == :als: use requested init directly (default multistart)
    • solver != :als: use BTDALSWarmStartInit(...) so first-order methods start from a good BTD-ALS warm point
  4. If the warm-start rel-error exceeds warm_rel_error_gate, return the warm BTD-ALS result directly
  5. Otherwise solve through _solve_model(...)
  6. If solver != :als, optionally polish with BTD-ALS (btd_als_polish_maxiter)
  7. Convert to BTDResult

Polish step usefulness (brief):

  • Usually helpful for a small final rel_error reduction after RGD converges near a good basin.
  • Most useful for quality-focused runs (benchmarks, final fits).
  • Can be skipped for speed-sensitive runs (btd_als_polish_maxiter=0) when small extra gains are not worth runtime.

BTD-specific initialization options:

  • :hosvd: sequential block initialization on residual
  • :hosvd_multistart: HOSVD subspace split candidates, optional screening ALS, keep lowest-cost candidate
  • :alswarm: short BTD-ALS warm-start wrapper around base initializer

BTD-ALS stabilization behavior:

  • Tracks per-iteration fit change (|rel_t - rel_{t-1}|)
  • Detects stagnation when fit change is tiny but rel_error remains high
  • Can restart from fresh multistart pool (max_stagnation_restarts)
  • Reports true final Riemannian gradient norm (grad_norm) instead of a placeholder

Tucker (tucker)

tucker(A, ranks; method=...) does not use _solve_model. It dispatches directly to decomposition routines:

  • :sthosvd
  • :hooi

Generic approx(...) routing

approx(manifolds, target; dispatch=:auto) routes by manifold family:

  • uniform Manifolds.Segre -> cpd(...)
  • uniform Manifolds.Tucker matching target shape/rank -> btd(...)
  • mixed or non-uniform family -> generic JoinModel(...) path -> ApproxResult

dispatch=:cpd, :btd, and :generic force behavior.

Result types and post-processing

  • CPDResult
  • BTDResult
  • TuckerResult
  • ApproxResult

Common utilities:

  • reconstruct(result)
  • rel_error(A, result)

File map

  • API entry points: src/api/approx.jl, src/api/cpd.jl, src/api/nncpd.jl, src/api/btd.jl
  • Routing helpers: src/dispatch/approx_routing.jl
  • Solver dispatch core: src/solvers/solve_dispatch.jl
  • BTD backend/init details: src/btd/model.jl, src/solvers/btd_als.jl
+Pipeline · TensorKitchen.jl
GitHub

TensorKitchen Pipeline

This document explains how public APIs route into models, solvers, and result converters.

Public entry points

  • CP Decomposition cpd(A, r; ...)
  • Nonnegative CP Decomposition nncpd(A, r; ...)
  • Block Term Decomposition btd(A, blocks, ranks; ...)
  • Tucker Decomposition tucker(A, ranks; method=...)
  • Join Decomposition approx(...)

Default behavior (quick reference)

  • cpd(A, r):
    • init = :alswarm
    • solver = :rgd
  • nncpd(A, r):
    • init = :alswarm
    • solver = :rgd
  • btd(A, blocks, ranks):
    • init = :alswarm
    • warm_steps = 200
    • warm_init = BTDHOSVDMultistartInit(candidates=64, screening_steps=10, block_maxiter=12)
    • warm_rel_error_gate = nothing (run manifold refinement by default; set e.g. 5e-2 to short-circuit on poor warm starts)
    • solver = :rgd
    • final BTD-ALS polish enabled by default for non-ALS solvers
    • max_stagnation_restarts = 1 (retry with stronger multistart when ALS fit-change stalls at high rel-error)
  • tucker(A, ranks):
    • method = :sthosvd
  • approx(model::JoinModel):
    • init = :random
    • solver = :rgd

Core execution architecture

Most optimization APIs share this core pattern:

  1. Build a model (JoinModel + backend)
  2. Call _solve_model(...)
  3. Convert to a public result struct

_solve_model lives in src/solvers/solve_dispatch.jl and is the common symbol-to-solver dispatch layer (:rgd, :rcg, :lbfgs, :als, :btd_tsd).

API flows

CPD (cpd, nncpd)

cpd(A, r; ...):

  1. Build JoinModel(A, r; geometry=...) with CPDBackend
  2. Normalize/validate options (solver, geometry, gradient_mode, normalization policy)
  3. Solve through _solve_model(...)
  4. Convert to CPDResult

Notes:

  • :als means CP-ALS.
  • Manifold solvers (:rgd, :rgd_fixed, :rcg, :lbfgs) share dispatch with other pipelines.
  • For solver != :als, init = :auto resolves to :alswarm, so CPD and NNCPD start from an ALS warm point before manifold refinement.
  • Generic approx(...) does not use CPD's ALS warm-start path unless it auto-routes to cpd(...).

BTD (btd)

btd(A, blocks, ranks; ...):

  1. Build a uniform Tucker family via TuckerJoin(...)
  2. Wrap as JoinModel with BTDBackend
  3. Choose effective initializer:
    • solver == :als: use requested init directly (default multistart)
    • solver != :als: use BTDALSWarmStartInit(...) so first-order methods start from a good BTD-ALS warm point
  4. If the warm-start rel-error exceeds warm_rel_error_gate, return the warm BTD-ALS result directly
  5. Otherwise solve through _solve_model(...)
  6. If solver != :als, optionally polish with BTD-ALS (btd_als_polish_maxiter)
  7. Convert to BTDResult

Polish step usefulness:

  • Usually helpful for a small final rel_error reduction after RGD converges near a good basin.
  • Most useful for quality-focused runs (benchmarks, final fits).
  • Can be skipped for speed-sensitive runs (btd_als_polish_maxiter=0) when small extra gains are not worth runtime.

BTD-specific initialization options:

  • :hosvd: sequential block initialization on residual
  • :hosvd_multistart: HOSVD subspace split candidates, optional screening ALS, keep lowest-cost candidate
  • :alswarm: short BTD-ALS warm-start wrapper around base initializer

BTD-ALS stabilization behavior:

  • Tracks per-iteration fit change (|rel_t - rel_{t-1}|)
  • Detects stagnation when fit change is tiny but rel_error remains high
  • Can restart from fresh multistart pool (max_stagnation_restarts)
  • Reports true final Riemannian gradient norm (grad_norm)

Tucker (tucker)

tucker(A, ranks; method=...) does not use _solve_model. It dispatches directly to decomposition routines:

  • :sthosvd
  • :hooi

Generic approx(...) routing

approx(manifolds, target; dispatch=:auto) routes by manifold family:

  • uniform Manifolds.Segre -> cpd(...)
  • uniform Manifolds.Tucker matching target shape/rank -> btd(...)
  • mixed or non-uniform family -> generic JoinModel(...) path -> ApproxResult

dispatch=:cpd, :btd, and :generic force behavior.

For the generic JoinModel path, approx(...) starts from init = :random by default and then runs the selected manifold solver. It does not run an ALS warm-start stage, because a general join component does not necessarily expose factor matrices or least-squares block updates.

Result types and post-processing

  • CPDResult
  • BTDResult
  • TuckerResult
  • ApproxResult

Common utilities:

  • reconstruct(result)
  • rel_error(A, result)

File map

  • API entry points: src/api/approx.jl, src/api/cpd.jl, src/api/nncpd.jl, src/api/btd.jl, src/api/tucker.jl
diff --git a/docs/btd/index.html b/docs/btd/index.html index 31fb97c..b214418 100644 --- a/docs/btd/index.html +++ b/docs/btd/index.html @@ -1,5 +1,5 @@ -BTD · Documentation
GitHub

Block Term Decomposition

A block term decomposition (BTD) with r blocks writes

\[\hat A = \sum_{i=1}^r A_i,\]

where each block $A_i$ is represented as a Tucker decomposition. At present, only homogeneous BTDs are supported, that is, all blocks must have the same multilinear rank.

To compute a block term decomposition of A with 10 blocks, each of multilinear rank (5, 4, 3), use

julia> r = 10
+BTD · TensorKitchen.jl
GitHub
Block Term Decomposition
A block term decomposition (BTD) with r blocks writes
\[\hat A = \sum_{i=1}^r A_i,\]
where each block $A_i$ is represented as a Tucker decomposition. At present, only homogeneous BTDs are supported, that is, all blocks must have the same multilinear rank.
To compute a block term decomposition of A with 10 blocks, each of multilinear rank (5, 4, 3), use
julia> r = 10
 julia> mlrank = (5, 4, 3)
 julia> btd_res = btd(A, r, mlrank)
 BTDResult{Float64}
@@ -12,6 +12,4 @@
 julia> res = btd(A, blocks, ranks; verbose = false)
 BTDResult{Float64}
   Blocks:       10
-  Rel. error:   0.2625821087015455
source
TensorKitchen.BTDResult — Type
BTDResult{T}

Result of block-term decomposition (btd); block components expose Tucker structure through accessors like core(blk), factors(blk), and blk.tensor.

  • solver_info: solver-specific diagnostics/metadata (NamedTuple). Typical keys include BTD-ALS restart diagnostics (total_iterations, stagnation_restarts, restart_rel_error_history) and BTD-TSD run settings (schedule, block_repeats, block_count, stepsize).
source
TensorKitchen.blocks — Function
blocks(r::BTDResult)

Return the Tucker block components of a block-term decomposition result.

source
TensorKitchen.reconstruct — Method
reconstruct(td::TuckerResult) reconstructs the tensor from Tucker decomposition
-
-A = S ×₁ U₁ ×₂ U₂ ⋯ ×_d U_d
source
+ Rel. error: 0.2625821087015455source
TensorKitchen.BTDResult — Type
BTDResult{T}

Result of block-term decomposition (btd); block components expose Tucker structure through accessors like core(blk), factors(blk), and blk.tensor.

  • solver_info: solver-specific diagnostics/metadata (NamedTuple). Typical keys include BTD-ALS restart diagnostics (total_iterations, stagnation_restarts, restart_rel_error_history) and BTD-TSD run settings (schedule, block_repeats, block_count, stepsize).
source
TensorKitchen.blocks — Function
blocks(r::BTDResult)

Return the Tucker block components of a block-term decomposition result.

source
TensorKitchen.reconstruct — Method
reconstruct(res::BTDResult)

Reconstruct the dense tensor represented by a block-term decomposition result by summing the reconstructed Tucker blocks.

source
diff --git a/docs/cpd/index.html b/docs/cpd/index.html index badcf94..9aee07a 100644 --- a/docs/cpd/index.html +++ b/docs/cpd/index.html @@ -1,5 +1,5 @@ -CPD · Documentation
GitHub

CPD

Here is how to approximate a tensor A by a CPD of rank r.

julia> using TensorKitchen
+CPD · TensorKitchen.jl
GitHub
CPD
Here is how to approximate a tensor A by a CPD of rank r.
julia> using TensorKitchen
 julia> A = randn(20, 15, 10)
 julia> r = 35
 julia> res = cpd(A, r)
@@ -14,13 +14,11 @@
   Order:        3
   Dimensions:   (20, 15, 10)
   Rank:         35
-  Rel. error:   0.4359141301703327
source
TensorKitchen.nncpd — Function
 nncpd(A, r; kwargs...)

Computes a nonnegative rank-r CP approximation of A in two steps: (1) the first step finds an initial point; (2) the second step refines the initial point. Returns a CPDResult. If r is omitted, uses the smallest tensor mode as a heuristic rank. cpd(A, r; nonnegative=true, ...) routes here and adopts the same effective defaults.

Options

The options are the same as for cpd.

Geometry guide:

  • geometry=:softplus_metric Default and usually the safest choice.
  • geometry=:squaring_metric Uses a regularized pullback-inspired geometry induced by the squaring chart.
  • geometry=:canonical Plain nonnegative CP coordinates without the pullback-style manifold geometry. This is the natural choice with solver=:als.

Example

julia> A = randn(20, 15, 10); r = 35;
+  Rel. error:   0.4359141301703327
source
TensorKitchen.nncpd — Function
 nncpd(A, r; kwargs...)

Computes a nonnegative rank-r CP approximation of A in two steps: (1) the first step finds an initial point; (2) the second step refines the initial point. Returns a CPDResult. If r is omitted, uses the smallest tensor mode as a heuristic rank. cpd(A, r; nonnegative=true, ...) routes here and adopts the same effective defaults.

Options

The options are the same as for cpd.

Geometry guide:

  • geometry=:softplus_metric Default and usually the safest choice.
  • geometry=:squaring_metric Uses a regularized pullback-inspired geometry induced by the squaring chart.
  • geometry=:canonical Plain nonnegative CP coordinates without the pullback-style manifold geometry. This is the natural choice with solver=:als.

Example

julia> A = randn(20, 15, 10); r = 35;
 julia> B = abs.(A)
 julia> nncpd(B, r)
 CPDResult{Float64}
   Order:        3
   Dimensions:   (20, 15, 10)
   Rank:         35
-  Rel. error:   0.3765605093526155
source
TensorKitchen.CPDResult — Type
CPDResult{T}

Result of a Canonical Polyadic Decomposition.

Stores the decoded CP representation together with solver diagnostics:

  • components: rank-one tensor components
  • weights: component weights
  • factors: factor matrices
  • cost: final objective function value at the returned solution
  • rel_error: final relative reconstruction error
  • grad_norm: norm of the final optimization gradient reported by the solver; for manifold solvers this is the Riemannian gradient norm
  • iterations: number of refinement iterations
  • converged: whether the solver reported convergence
  • solver: solver optimization method used to produce the result
  • solver_info: solver-specific diagnostics/metadata (NamedTuple). Typical keys include:
    • initial_stepsize_eff (RGD), memory_size (LBFGS), cautious_update (LBFGS),
    • initial_scale, linesearch, has_preconditioner (LBFGS), and nncp_pullback_eps (NNCP).
source
TensorKitchen.factors — Method
factors(res::CPDResult)

Return the CP factor matrices of res as a vector [U₁, U₂, ..., U_N], where each U_m has size size(A, m) × rank.

source
TensorKitchen.reconstruct — Method
reconstruct(td::TuckerResult) reconstructs the tensor from Tucker decomposition
-
-A = S ×₁ U₁ ×₂ U₂ ⋯ ×_d U_d
source
+ Rel. error: 0.3765605093526155source
TensorKitchen.CPDResult — Type
CPDResult{T}

Result of a Canonical Polyadic Decomposition.

Stores the decoded CP representation together with solver diagnostics:

  • components: rank-one tensor components
  • weights: component weights
  • factors: factor matrices
  • cost: final objective function value at the returned solution
  • rel_error: final relative reconstruction error
  • grad_norm: norm of the final optimization gradient reported by the solver; for manifold solvers this is the Riemannian gradient norm
  • iterations: number of refinement iterations
  • converged: whether the solver reported convergence
  • solver: solver optimization method used to produce the result
  • solver_info: solver-specific diagnostics/metadata (NamedTuple). Typical keys include:
    • initial_stepsize_eff (RGD), memory_size (LBFGS), cautious_update (LBFGS),
    • initial_scale, linesearch, has_preconditioner (LBFGS), and nncp_pullback_eps (NNCP).
source
TensorKitchen.weights — Method
weights(r::CPDResult)

Return the CP component weights stored in a CPD result.

source
TensorKitchen.factors — Method
factors(res::CPDResult)

Return the CP factor matrices of res as a vector [U₁, U₂, ..., U_N], where each U_m has size size(A, m) × rank.

source
TensorKitchen.reconstruct — Method
reconstruct(res::CPDResult)

Reconstruct the dense tensor represented by a CP decomposition result.

For a rank-R CPD result, this returns sum(weights(res)[k] * u_1k ⊗ ... ⊗ u_Nk for k = 1:R).

source
diff --git a/docs/index.html b/docs/index.html index fbd74d5..6a69566 100644 --- a/docs/index.html +++ b/docs/index.html @@ -1,2 +1,2 @@ -Home · Documentation
GitHub

TensorKitchen.jl Documentation

TensorKitchen.jl is a Julia package for tensor decompositions.

Notes

The package is currently at a pre-alpha stage.

The implementation is based on combining algebraic algorithms like ALS (see, e.g., the textbook by Kolda and Ballard) and Riemannian optimization from Manopt.jl.

What currently works is

  • Canonical Polyadic Decomposition (CPD)
  • Tucker Decomposition
  • Nonnegative Canonical Polyadic Decomposition (NNCPD)
  • Block Term Decomposition (BTD)
  • Join Decompositions

The next updates will include

  • Handling of swamps/plateaus in the optimization step
  • Documentation
  • Improved User Interface
  • GPU Support
  • LL1 Decomposition (3-way specialized BTD)
  • Symmetric CP / Waring Decomposition
  • Partially Symmetric CP
  • Tensor Trains
+Home · TensorKitchen.jl
GitHub

TensorKitchen.jl: tensor decompositions in Julia

<img src="logotexttransparent.png" width="450px">

TensorKitchen.jl is a Julia package for tensor decompositions.

Notes

The package is currently an early version and will be updated frequently in the near future.

The implementation is based on combining algebraic algorithms like ALS (see, e.g., the textbook by Kolda and Ballard) and Riemannian optimization from Manopt.jl.

What currently works is

  • Canonical Polyadic Decomposition (CPD)
  • Tucker Decomposition
  • Nonnegative Canonical Polyadic Decomposition (NNCPD)
  • Block Term Decomposition (BTD)
  • Join Decompositions

The next updates will include

  • Handling of swamps/plateaus in the optimization step
  • Documentation
  • Improved User Interface
  • GPU Support
  • LL1 Decomposition (3-way specialized BTD)
  • Symmetric CP / Waring Decomposition
  • Partially Symmetric CP
  • Tensor Trains
diff --git a/docs/join/index.html b/docs/join/index.html index 80a16d8..cded9ca 100644 --- a/docs/join/index.html +++ b/docs/join/index.html @@ -1,5 +1,5 @@ -Join · Documentation
GitHub

Join Decomposition

A Join Decomposition of a vector $x\in\mathbb R^N$ is a decomposition of the form $x = x_1+\cdots+x_r$, where $x_i\in M_i$ and $M_i\subset \mathbb R^N$ is a given embedded manifold.

For instance, we can approximate a point $p = (1.2, 0.4)\in\mathbb R^2$ by $x=x_1+x_2$, where $x_1,x_2\in S^1$ are points on the circle:

julia> using Manifolds
+Join · TensorKitchen.jl
GitHub
Join Decomposition
A Join Decomposition of a vector $x\in\mathbb R^N$ is a decomposition of the form $x = x_1+\cdots+x_r$, where $x_i\in M_i$ and $M_i\subset \mathbb R^N$ is a given embedded manifold. 
For instance, we can approximate a point $p = (1.2, 0.4)\in\mathbb R^2$ by $x=x_1+x_2$, where $x_1,x_2\in S^1$ are points on the circle:
julia> using Manifolds
 julia> p = [1.2, 0.4]
 julia> S = Sphere(1)
 julia> join_res = approx((S, S), p)
@@ -37,5 +37,4 @@
 approx(Manifolds.Segre((2, 3)), 2, target; verbose = false)
 # returns an CPDResult
  • approx(base, target; kwargs...) : Builds a single-component generic join and route to the generic join solver.   Example:
target = [1.2, 0.4, -0.3]   
 approx(Manifolds.Sphere(2), target; verbose = false)
-# returns an ApproxResult
  • By default, approx auto-routes by manifold family:
    • uniform Manifolds.Segre summands calls cpd(...)
    • uniform Manifolds.Tucker summands calls btd(...)
    • otherwise calls JoinModel(...) and returns a ApproxResult
Return Types
  • Depending on the manifold family, approx(...) may return:
    • ApproxResult for the generic join path
    • CPDResult when auto-routed to cpd(...)
    • BTDResult when auto-routed to btd(...)
Main Options
For the generic join path:
  • init = :random: Sets the algorithm to find the initial point.
  • solver = :rgd: Sets the algorithm for refinement. Possible options are:
    • rgd (default): Riemannian gradient descent
    • rgd_fixed: Riemannian gradient descent with fixed step size
    • rcg: Riemannian conjugate gradient
    • lbfgs: Limited-memory quasi-Newton
##Notes##
  • :als is not a solver option for approx(...). However, if approx(...) auto-routes to cpd(...) or btd(...), then those specialized pipelines may support ALS separately.
  • warm_steps and warm_init are not part of the generic approx(...) path. Generic joins start from the selected initializer and then use manifold solvers for refinement.
  • For generic mixed joins, use manifold solvers such as :rgd, :rcg, or :lbfgs.
source
approx(M::ProductManifold, target; dispatch=:auto, kwargs...)
Use the factors of a product manifold as join components and route to CPD, BTD, or the generic join solver according to dispatch.
source
approx(base::Manifolds.Segre, r, target; dispatch=:auto, kwargs...)
Build a rank-r Segre join and route to the CPD pipeline unless generic dispatch is explicitly requested.
source
approx(base::Manifolds.Tucker, r, target; dispatch=:auto, kwargs...)
Build a r-block Tucker join and route to the BTD pipeline unless generic dispatch is explicitly requested.
source
approx(base::AbstractManifold, r, target; dispatch=:auto, kwargs...)
Fallback rank-r join constructor for non-specialized manifolds. Forced CPD or BTD dispatch is rejected because the base manifold family is not known.
source
approx(base::AbstractManifold, target; dispatch=:auto, kwargs...)
Single-component generic approximation fallback. Use this when no CPD/BTD family-specific routing is intended.
source
TensorKitchen.ApproxResult — Type
ApproxResult{T}

Generic result of approx(manifolds, target) (generic join decomposition).

  • point: final point on the join manifold
  • components: extracted component descriptions
  • cost: the value of the optimization objective at the final returned point, that is typically the least-squares objective.
  • rel_error: relative error of the decomposition which is the ratio of the cost to the target norm. is scale-normalized and easier to compare across problems
  • grad_norm: norm of the final optimization gradient reported by the solver; for manifold solvers this is typically the Riemannian gradient norm
  • iterations: number of iterations used
  • converged: whether the decomposition converged
  • solver: solver used for the decomposition
  • solver_info: solver-specific diagnostics/metadata (NamedTuple)
source
TensorKitchen.join_product — Function
join_product(base, r) constructs a ProductManifold

Construct rank-r join parameter manifold as a plain ProductManifold. For Manifolds.Segre, uses flattened (Euclidean(1), Sphere, ..., Sphere) factors per component to preserve current CP parameter layout.

source
TensorKitchen.SegreProduct — Function
SegreProduct(dims, r)

Product manifold Manifolds.Segre(dims) × ... × Manifolds.Segre(dims) (r factors). Each component point uses the Manifolds.Segre layout [[λ], x₁, …, x_d].

source
+# returns an ApproxResult

Return Types

Main Options

For the generic join path:

##Notes##

source
approx(M::ProductManifold, target; dispatch=:auto, kwargs...)

Use the factors of a product manifold as join components and route to CPD, BTD, or the generic join solver according to dispatch.

source
approx(base::Manifolds.Segre, r, target; dispatch=:auto, kwargs...)

Build a rank-r Segre join and route to the CPD pipeline unless generic dispatch is explicitly requested.

source
approx(base::Manifolds.Tucker, r, target; dispatch=:auto, kwargs...)

Build a r-block Tucker join and route to the BTD pipeline unless generic dispatch is explicitly requested.

source
approx(base::AbstractManifold, r, target; dispatch=:auto, kwargs...)

Fallback rank-r join constructor for non-specialized manifolds. Forced CPD or BTD dispatch is rejected because the base manifold family is not known.

source
approx(base::AbstractManifold, target; dispatch=:auto, kwargs...)

Single-component generic approximation fallback. Use this when no CPD/BTD family-specific routing is intended.

source
TensorKitchen.ApproxResult — Type
ApproxResult{T}

Generic result of approx(manifolds, target) (generic join decomposition).

  • point: final point on the join manifold
  • components: extracted component descriptions
  • cost: the value of the optimization objective at the final returned point, that is typically the least-squares objective.
  • rel_error: relative error of the decomposition which is the ratio of the cost to the target norm. is scale-normalized and easier to compare across problems
  • grad_norm: norm of the final optimization gradient reported by the solver; for manifold solvers this is typically the Riemannian gradient norm
  • iterations: number of iterations used
  • converged: whether the decomposition converged
  • solver: solver used for the decomposition
  • solver_info: solver-specific diagnostics/metadata (NamedTuple)
source
TensorKitchen.reconstruct — Method
reconstruct(res::ApproxResult)

Reconstruct the dense ambient object represented by a generic join approximation result by summing its component tensors.

source
TensorKitchen.join_product — Function
join_product(base, r)

Decides how to expand base into r components:

  • Manifolds.Segre → flattened (Euclidean(1), Sphere, ...) × r (one λ + spheres per rank-1).
  • Manifolds.Tucker / ProductManifold → repeat each factor r times.
  • Generic manifold → ProductManifold(base, base, ..., base).
source
TensorKitchen.SegreProduct — Function
SegreProduct(dims, r)
  • Product manifold Manifolds.Segre(dims) × ... × Manifolds.Segre(dims) (r factors).
  • Each component point uses the Manifolds.Segre layout [[λ], x₁, …, x_d].
  • SegreProduct is used to build a join model for CPD.
source
diff --git a/docs/logo_text_transparent.png b/docs/logo_text_transparent.png index 5e114b1..f12234f 100644 Binary files a/docs/logo_text_transparent.png and b/docs/logo_text_transparent.png differ diff --git a/docs/make.jl b/docs/make.jl index c8e1a33..5a3c770 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,7 +1,7 @@ using Documenter, TensorKitchen makedocs( - sitename = "Documentation", + sitename = "TensorKitchen.jl", warnonly = true, pages = [ "Home" => "index.md", diff --git a/docs/references/index.html b/docs/references/index.html index 46f79a8..0b4dc7c 100644 --- a/docs/references/index.html +++ b/docs/references/index.html @@ -1,2 +1,2 @@ -References · Documentation
GitHub

References

General Tensor Decomposition

  • Tensor decompositions (CP, Tucker): T. G. Kolda and B. W. Bader, "Tensor decompositions and applications," SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.

Tucker Methods

  • HOSVD: L. De Lathauwer, B. De Moor, and J. Vandewalle, "A multilinear singular value decomposition," SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1253–1278, 2000.
  • ST-HOSVD: N. Vannieuwenhoven, R. Vandebril, K. Meerbergen, "A new truncation strategy for the higher-order singular value decomposition," SIAM J. Sci. Comput., vol. 34, no. 2, pp. A1027–A1052, 2012.
  • HOOI: L. De Lathauwer, B. De Moor, and J. Vandewalle, "On the best rank-1 and rank-(R1,R2,...,R_N) approximation of higher-order tensors," SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1324–1342, 2000.

Block and Structured Models (BTD / LL1)

  • Block-term decomposition (BTD): L. De Lathauwer, "Decompositions of a higher-order tensor in block terms—Part I: Lemmas for partitioned matrices," SIAM J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1022–1032, 2008.
  • L. De Lathauwer, "Decompositions of a higher-order tensor in block terms—Part II: Definitions and uniqueness," SIAM J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1033–1066, 2008.
  • BTD-ALS: L. De Lathauwer and D. Nion, "Decompositions of a higher-order tensor in block terms—Part III: Alternating least squares algorithms," SIAM Journal on Matrix Analysis and Applications, vol. 30, no. 3, pp. 1067–1083, 2008. PDF.

Join decompositions

  • Conditioning of join decompositions: P. Breiding and N. Vannieuwenhoven, "The condition number of join decompositions," SIAM Journal on Matrix Analysis and Applications, vol. 39, no. 1, pp. 287–309, 2018. arXiv:1611.08117 (PDF).

Riemannian Optimization and Julia Ecosystem

  • Riemannian trust-region / Gauss–Newton for canonical rank (CP) approximation: P. Breiding and N. Vannieuwenhoven, "A Riemannian Trust Region Method for the Canonical Tensor Rank Approximation Problem," SIAM Journal on Optimization, vol. 28, no. 3, pp. 2435–2465, 2018. arXiv:1709.00033 (PDF).
  • Riemannian optimization: P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008.
  • Julia manifold optimization ecosystem: R. Bergmann et al., ManifoldsBase.jl, Manifolds.jl, and Manopt.jl.
+References · TensorKitchen.jl
GitHub

References

General Tensor Decomposition

  • Tensor decompositions (CP, Tucker): T. G. Kolda and B. W. Bader, "Tensor decompositions and applications," SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.

Tucker Methods

  • HOSVD: L. De Lathauwer, B. De Moor, and J. Vandewalle, "A multilinear singular value decomposition," SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1253–1278, 2000.
  • ST-HOSVD: N. Vannieuwenhoven, R. Vandebril, K. Meerbergen, "A new truncation strategy for the higher-order singular value decomposition," SIAM J. Sci. Comput., vol. 34, no. 2, pp. A1027–A1052, 2012.
  • HOOI: L. De Lathauwer, B. De Moor, and J. Vandewalle, "On the best rank-1 and rank-(R1,R2,...,R_N) approximation of higher-order tensors," SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1324–1342, 2000.

Block and Structured Models (BTD / LL1)

  • Block-term decomposition (BTD): L. De Lathauwer, "Decompositions of a higher-order tensor in block terms—Part I: Lemmas for partitioned matrices," SIAM J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1022–1032, 2008.
  • L. De Lathauwer, "Decompositions of a higher-order tensor in block terms—Part II: Definitions and uniqueness," SIAM J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1033–1066, 2008.
  • BTD-ALS: L. De Lathauwer and D. Nion, "Decompositions of a higher-order tensor in block terms—Part III: Alternating least squares algorithms," SIAM Journal on Matrix Analysis and Applications, vol. 30, no. 3, pp. 1067–1083, 2008. PDF.

Join decompositions

  • Conditioning of join decompositions: P. Breiding and N. Vannieuwenhoven, "The condition number of join decompositions," SIAM Journal on Matrix Analysis and Applications, vol. 39, no. 1, pp. 287–309, 2018. arXiv:1611.08117 (PDF).

Riemannian Optimization and Julia Ecosystem

  • Riemannian trust-region / Gauss–Newton for canonical rank (CP) approximation: P. Breiding and N. Vannieuwenhoven, "A Riemannian Trust Region Method for the Canonical Tensor Rank Approximation Problem," SIAM Journal on Optimization, vol. 28, no. 3, pp. 2435–2465, 2018. arXiv:1709.00033 (PDF).
  • Riemannian optimization: P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008.
  • Julia manifold optimization ecosystem: R. Bergmann et al., ManifoldsBase.jl, Manifolds.jl, and Manopt.jl.
diff --git a/docs/search_index.js b/docs/search_index.js index 4d9e21c..efd178b 100644 --- a/docs/search_index.js +++ b/docs/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"category":"section","location":"references/#References","page":"References","text":"","title":"References"},{"category":"section","location":"references/#General-Tensor-Decomposition","page":"References","text":"Tensor decompositions (CP, Tucker): T. G. Kolda and B. W. Bader, \"Tensor decompositions and applications,\" SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.","title":"General Tensor Decomposition"},{"category":"section","location":"references/#Tucker-Methods","page":"References","text":"HOSVD: L. De Lathauwer, B. De Moor, and J. Vandewalle, \"A multilinear singular value decomposition,\" SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1253–1278, 2000.\nST-HOSVD: N. Vannieuwenhoven, R. Vandebril, K. Meerbergen, \"A new truncation strategy for the higher-order singular value decomposition,\" SIAM J. Sci. Comput., vol. 34, no. 2, pp. A1027–A1052, 2012.\nHOOI: L. De Lathauwer, B. De Moor, and J. Vandewalle, \"On the best rank-1 and rank-(R1,R2,...,R_N) approximation of higher-order tensors,\" SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1324–1342, 2000.","title":"Tucker Methods"},{"category":"section","location":"references/#Block-and-Structured-Models-(BTD-/-LL1)","page":"References","text":"Block-term decomposition (BTD): L. De Lathauwer, \"Decompositions of a higher-order tensor in block terms—Part I: Lemmas for partitioned matrices,\" SIAM J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1022–1032, 2008.\nL. De Lathauwer, \"Decompositions of a higher-order tensor in block terms—Part II: Definitions and uniqueness,\" SIAM J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1033–1066, 2008.\nBTD-ALS: L. De Lathauwer and D. Nion, \"Decompositions of a higher-order tensor in block terms—Part III: Alternating least squares algorithms,\" SIAM Journal on Matrix Analysis and Applications, vol. 30, no. 3, pp. 1067–1083, 2008. PDF.","title":"Block and Structured Models (BTD / LL1)"},{"category":"section","location":"references/#Join-decompositions","page":"References","text":"Conditioning of join decompositions: P. Breiding and N. Vannieuwenhoven, \"The condition number of join decompositions,\" SIAM Journal on Matrix Analysis and Applications, vol. 39, no. 1, pp. 287–309, 2018. arXiv:1611.08117 (PDF).","title":"Join decompositions"},{"category":"section","location":"references/#Riemannian-Optimization-and-Julia-Ecosystem","page":"References","text":"Riemannian trust-region / Gauss–Newton for canonical rank (CP) approximation: P. Breiding and N. Vannieuwenhoven, \"A Riemannian Trust Region Method for the Canonical Tensor Rank Approximation Problem,\" SIAM Journal on Optimization, vol. 28, no. 3, pp. 2435–2465, 2018. arXiv:1709.00033 (PDF).\nRiemannian optimization: P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008.\nJulia manifold optimization ecosystem: R. Bergmann et al., ManifoldsBase.jl, Manifolds.jl, and Manopt.jl.","title":"Riemannian Optimization and Julia Ecosystem"},{"category":"section","location":"PIPELINE/#TensorKitchen-Pipeline","page":"Pipeline","text":"This document explains how public APIs route into models, solvers, and result converters.","title":"TensorKitchen Pipeline"},{"category":"section","location":"PIPELINE/#Public-entry-points","page":"Pipeline","text":"cpd(A, r; ...) -> CPDResult\nnncpd(A, r; ...) -> CPDResult\nbtd(A, blocks, ranks; ...) -> BTDResult\ntucker(A, ranks; method=...) -> TuckerResult\napprox(...) -> ApproxResult or auto-routed CPDResult/BTDResult","title":"Public entry points"},{"category":"section","location":"PIPELINE/#Default-behavior-(quick-reference)","page":"Pipeline","text":"cpd(A, r):\ninit = :alswarm\nsolver = :rgd\nnncpd(A, r):\ninit = :alswarm\nsolver = :rgd\nbtd(A, blocks, ranks):\ninit = :alswarm\nwarm_steps = 200\nwarm_init = BTDHOSVDMultistartInit(candidates=64, screening_steps=10, block_maxiter=12)\nwarm_rel_error_gate = nothing (run manifold refinement by default; set e.g. 5e-2 to short-circuit on poor warm starts)\nsolver = :rgd\nfinal BTD-ALS polish enabled by default for non-ALS solvers\nmax_stagnation_restarts = 1 (retry with stronger multistart when ALS fit-change stalls at high rel-error)\ntucker(A, ranks):\nmethod = :sthosvd\napprox(model::JoinModel):\ninit = :alswarm\nwarm_steps = 500\nsolver = :rgd","title":"Default behavior (quick reference)"},{"category":"section","location":"PIPELINE/#Core-execution-architecture","page":"Pipeline","text":"Most optimization APIs share this core pattern:\n\nBuild a model (JoinModel + backend)\nCall _solve_model(...)\nConvert to a public result struct\n\n_solve_model lives in src/solvers/solve_dispatch.jl and is the common symbol-to-solver dispatch layer (:rgd, :rcg, :lbfgs, :als, :btd_tsd).","title":"Core execution architecture"},{"category":"section","location":"PIPELINE/#API-flows","page":"Pipeline","text":"","title":"API flows"},{"category":"section","location":"PIPELINE/#CPD-(cpd,-nncpd)","page":"Pipeline","text":"cpd(A, r; ...):\n\nBuild JoinModel(A, r; geometry=...) with CPDBackend\nNormalize/validate options (solver, geometry, gradient_mode, normalization policy)\nSolve through _solve_model(...)\nOptionally run nonnegative ALS polishing (for selected nonnegative paths)\nConvert to CPDResult\n\nNotes:\n\n:als means CP-ALS.\nManifold solvers (:rgd, :rgd_fixed, :rcg, :lbfgs) share dispatch with other pipelines.","title":"CPD (cpd, nncpd)"},{"category":"section","location":"PIPELINE/#BTD-(btd)","page":"Pipeline","text":"btd(A, blocks, ranks; ...):\n\nBuild a uniform Tucker family via TuckerJoin(...)\nWrap as JoinModel with BTDBackend\nChoose effective initializer:\nsolver == :als: use requested init directly (default multistart)\nsolver != :als: use BTDALSWarmStartInit(...) so first-order methods start from a good BTD-ALS warm point\nIf the warm-start rel-error exceeds warm_rel_error_gate, return the warm BTD-ALS result directly\nOtherwise solve through _solve_model(...)\nIf solver != :als, optionally polish with BTD-ALS (btd_als_polish_maxiter)\nConvert to BTDResult\n\nPolish step usefulness (brief):\n\nUsually helpful for a small final rel_error reduction after RGD converges near a good basin.\nMost useful for quality-focused runs (benchmarks, final fits).\nCan be skipped for speed-sensitive runs (btd_als_polish_maxiter=0) when small extra gains are not worth runtime.\n\nBTD-specific initialization options:\n\n:hosvd: sequential block initialization on residual\n:hosvd_multistart: HOSVD subspace split candidates, optional screening ALS, keep lowest-cost candidate\n:alswarm: short BTD-ALS warm-start wrapper around base initializer\n\nBTD-ALS stabilization behavior:\n\nTracks per-iteration fit change (|rel_t - rel_{t-1}|)\nDetects stagnation when fit change is tiny but rel_error remains high\nCan restart from fresh multistart pool (max_stagnation_restarts)\nReports true final Riemannian gradient norm (grad_norm) instead of a placeholder","title":"BTD (btd)"},{"category":"section","location":"PIPELINE/#Tucker-(tucker)","page":"Pipeline","text":"tucker(A, ranks; method=...) does not use _solve_model. It dispatches directly to decomposition routines:\n\n:sthosvd\n:hooi","title":"Tucker (tucker)"},{"category":"section","location":"PIPELINE/#Generic-approx(...)-routing","page":"Pipeline","text":"approx(manifolds, target; dispatch=:auto) routes by manifold family:\n\nuniform Manifolds.Segre -> cpd(...)\nuniform Manifolds.Tucker matching target shape/rank -> btd(...)\nmixed or non-uniform family -> generic JoinModel(...) path -> ApproxResult\n\ndispatch=:cpd, :btd, and :generic force behavior.","title":"Generic approx(...) routing"},{"category":"section","location":"PIPELINE/#Result-types-and-post-processing","page":"Pipeline","text":"CPDResult\nBTDResult\nTuckerResult\nApproxResult\n\nCommon utilities:\n\nreconstruct(result)\nrel_error(A, result)","title":"Result types and post-processing"},{"category":"section","location":"PIPELINE/#File-map","page":"Pipeline","text":"API entry points: src/api/approx.jl, src/api/cpd.jl, src/api/nncpd.jl, src/api/btd.jl\nRouting helpers: src/dispatch/approx_routing.jl\nSolver dispatch core: src/solvers/solve_dispatch.jl\nBTD backend/init details: src/btd/model.jl, src/solvers/btd_als.jl","title":"File map"},{"category":"section","location":"join/#Join-Decomposition","page":"Join","text":"A Join Decomposition of a vector xinmathbb R^N is a decomposition of the form x = x_1+cdots+x_r, where x_iin M_i and M_isubset mathbb R^N is a given embedded manifold. \n\nFor instance, we can approximate a point p = (12 04)inmathbb R^2 by x=x_1+x_2, where x_1x_2in S^1 are points on the circle:\n\njulia> using Manifolds\njulia> p = [1.2, 0.4]\njulia> S = Sphere(1)\njulia> join_res = approx((S, S), p)\nApproxResult{Float64}\n Components: 2\n Rel. error: 0.0007114699550529457\n\nWe access the decomposition as follows.\n\ncomponents(join_res)\nreconstruct(join_res)\n\n
","title":"Join Decomposition"},{"category":"section","location":"join/#Generic-Join-Approximation","page":"Join","text":"approx(...) is the main frontend for join decomposition. It works in two stages:\n\nbuild an initial point\nrefine it with the selected solver","title":"Generic Join Approximation"},{"category":"section","location":"join/#Supported-Forms","page":"Join","text":"approx(model; kwargs...) Solve an already constructed JoinModel.\napprox(manifolds, target; kwargs...) Build a join from a tuple/vector of component manifolds.\napprox(M::ProductManifold, target; kwargs...) Use the product-manifold factors as join components.\napprox(base, r, target; kwargs...) Repeat a single base manifold r times to build a join.\napprox(base, target; kwargs...) Build a single-component join.\n\nExamples:\n\njulia> target = [1.2, 0.4, -0.3]\njulia> model = JoinModel(Manifolds.Sphere(2), target)\njulia> approx(model; maxiter = 50, verbose = false)\nApproxResult{Float64}\n Components: 1\n Rel. error: 0.23076923076923075\n\njulia> target = randn(2, 3)\njulia> approx((Manifolds.Segre((2, 3)), Manifolds.Segre((2, 3))), target; verbose = false)\nCPDResult{Float64}\n\njulia> target = [1.2, 0.4, -0.3]\njulia> approx(Manifolds.Sphere(2), 2, target; verbose = false)\nApproxResult{Float64}","title":"Supported Forms"},{"category":"section","location":"join/#Routing","page":"Join","text":"With dispatch = :auto:\n\nuniform Manifolds.Segre components route to cpd(...)\nuniform Manifolds.Tucker components route to btd(...)\notherwise the generic JoinModel(...) path is used and ApproxResult is returned\n\nYou can also force the route explicitly:\n\ndispatch = :generic\ndispatch = :cpd\ndispatch = :btd\n\nForced routing is validated:\n\ndispatch = :cpd requires all components to be Manifolds.Segre with identical factor_dims\ndispatch = :btd requires all components to be Manifolds.Tucker with identical factor_dims and compatible multilinear rank","title":"Routing"},{"category":"section","location":"join/#Generic-Join-Options","page":"Join","text":"For the generic join path:\n\ninit = :random Default initializer. Built-in initializer support depends on the component manifolds.\nsolver = :rgd Supported generic-join solver options are:\n:rgd\n:rgd_fixed\n:rcg\n:lbfgs\n\nOther common options:\n\np0 = nothing\nmaxiter = 500\nstepsize = 1.0\ntol = 1e-6\ngradient_mode = :riemannian\nverbose = true\nvector_transport_method = nothing\n\nBuilt-in init support by component family:\n\nSphere: :random, :deterministic, :target\nSegre: :random, :deterministic\nTucker: :random, :tucker, :tucker_diag, :sthosvd\nother manifolds: :random","title":"Generic Join Options"},{"category":"section","location":"join/#Notes","page":"Join","text":"Generic joins require every component manifold to embed into the same flattened ambient length as target.\nsolver = :als is not available for a truly generic JoinModel. ALS may still be available when approx(...) auto-routes to cpd(...) or btd(...).","title":"Notes"},{"category":"section","location":"join/#Join-Decomposition-Docs","page":"Join","text":"","title":"Join Decomposition Docs"},{"category":"function","location":"join/#TensorKitchen.approx","page":"Join","text":"Generic Join Approximation\n\napprox(...) is the main frontend for join decomposition, which works in two stages:\n\nbuild an initial point\nrefine it with the selected solver\n\nSupported Forms\n\napprox(model; kwargs...) : It is for an already constructed join model (JoinModel(...)) and routes to the generic join solver. Model can be a JoinModel of a tuple of manifolds, a ProductManifold, or a single manifold.\nmodel means an already constructed JoinModel.\nIt fully fixes the decomposition structure and target.\napprox(model; ...) just solves that model.\nExample: If you want to approximate a point on the sphere, you can build a JoinModel and then use approx to refine it.\n\ntarget = [1.2, 0.4, -0.3]\nmodel = JoinModel(Manifolds.Sphere(2), target)\napprox(model; maxiter = 100, verbose = false)\n# returns an ApproxResult\n\napprox(manifolds, target; kwargs...) : Builds a generic join model from a tuple of existing manifolds and routes to according to dispatch. Example:\n\ntarget = randn(2, 3)\napprox((Manifolds.Segre((2, 3)), Manifolds.Segre((2, 3))), target; verbose = false)\n# This builds a join with 2 copies of Manifolds.Segre((2, 3))\".\n\napprox(M::ProductManifold, target; kwargs...) : Uses the factors of a product manifold and route to CPD, BTD, or the generic join solver according to dispatch. \nM::ProductManifold means that you already have a product manifold whose factors are the join components.\napprox(M, target; ...) uses those factors directly.\nIt is more explicit than base, because the components are already listed.\nExample:\n\n# Example 1\ntarget = randn(2, 3)\nM = ProductManifold(Manifolds.Segre((2, 3)), Manifolds.Segre((2, 3)))\napprox(M, target; verbose = false)\n# returns an CPDResult\n\n# Example 2\ntarget = randn(4, 3, 2)\nM = ProductManifold(\n Manifolds.Tucker((4, 3, 2), (2, 2, 2)),\n Manifolds.Tucker((4, 3, 2), (2, 2, 2)),\n)\napprox(M, target; verbose = false)\n# returns an BTDResult\n\napprox(base, r, target; kwargs...) : builds a rank-r Segre join and routes to CPD when base isa Manifolds.Segre and BTD when base isa Manifolds.Tucker unless generic\n\ndispatch is explicitly requested. - base means one manifold template, not yet a full join. - approx(base, r, target; ...) repeats that same manifold r times to build a join. - approx(base, target; ...) builds a one-component join.\n\ntarget = randn(2, 3)\napprox(Manifolds.Segre((2, 3)), 2, target; verbose = false)\n# returns an CPDResult\n\napprox(base, target; kwargs...) : Builds a single-component generic join and route to the generic join solver. Example:\n\ntarget = [1.2, 0.4, -0.3] \napprox(Manifolds.Sphere(2), target; verbose = false)\n# returns an ApproxResult\n\nBy default, approx auto-routes by manifold family:\nuniform Manifolds.Segre summands calls cpd(...)\nuniform Manifolds.Tucker summands calls btd(...)\notherwise calls JoinModel(...) and returns a ApproxResult\n\nReturn Types\n\nDepending on the manifold family, approx(...) may return:\nApproxResult for the generic join path\nCPDResult when auto-routed to cpd(...)\nBTDResult when auto-routed to btd(...)\n\nMain Options\n\nFor the generic join path:\n\ninit = :random: Sets the algorithm to find the initial point.\nsolver = :rgd: Sets the algorithm for refinement. Possible options are:\nrgd (default): Riemannian gradient descent\nrgd_fixed: Riemannian gradient descent with fixed step size\nrcg: Riemannian conjugate gradient\nlbfgs: Limited-memory quasi-Newton\n\n##Notes##\n\n:als is not a solver option for approx(...). However, if approx(...) auto-routes to cpd(...) or btd(...), then those specialized pipelines may support ALS separately.\nwarm_steps and warm_init are not part of the generic approx(...) path. Generic joins start from the selected initializer and then use manifold solvers for refinement.\nFor generic mixed joins, use manifold solvers such as :rgd, :rcg, or :lbfgs.\n\n\n\n\n\napprox(M::ProductManifold, target; dispatch=:auto, kwargs...)\n\nUse the factors of a product manifold as join components and route to CPD, BTD, or the generic join solver according to dispatch.\n\n\n\n\n\napprox(base::Manifolds.Segre, r, target; dispatch=:auto, kwargs...)\n\nBuild a rank-r Segre join and route to the CPD pipeline unless generic dispatch is explicitly requested.\n\n\n\n\n\napprox(base::Manifolds.Tucker, r, target; dispatch=:auto, kwargs...)\n\nBuild a r-block Tucker join and route to the BTD pipeline unless generic dispatch is explicitly requested.\n\n\n\n\n\napprox(base::AbstractManifold, r, target; dispatch=:auto, kwargs...)\n\nFallback rank-r join constructor for non-specialized manifolds. Forced CPD or BTD dispatch is rejected because the base manifold family is not known.\n\n\n\n\n\napprox(base::AbstractManifold, target; dispatch=:auto, kwargs...)\n\nSingle-component generic approximation fallback. Use this when no CPD/BTD family-specific routing is intended.\n\n\n\n\n\n","title":"TensorKitchen.approx"},{"category":"type","location":"join/#TensorKitchen.ApproxResult","page":"Join","text":"ApproxResult{T}\n\nGeneric result of approx(manifolds, target) (generic join decomposition).\n\npoint: final point on the join manifold\ncomponents: extracted component descriptions\ncost: the value of the optimization objective at the final returned point, that is typically the least-squares objective.\nrel_error: relative error of the decomposition which is the ratio of the cost to the target norm. is scale-normalized and easier to compare across problems\ngrad_norm: norm of the final optimization gradient reported by the solver; for manifold solvers this is typically the Riemannian gradient norm\niterations: number of iterations used\nconverged: whether the decomposition converged\nsolver: solver used for the decomposition\nsolver_info: solver-specific diagnostics/metadata (NamedTuple)\n\n\n\n\n\n","title":"TensorKitchen.ApproxResult"},{"category":"method","location":"join/#TensorKitchen.reconstruct-Tuple{ApproxResult}","page":"Join","text":"reconstruct(res::ApproxResult)\nreconstruct(res::BTDResult)\n\n\n\n\n\n","title":"TensorKitchen.reconstruct"},{"category":"function","location":"join/#TensorKitchen.join_product","page":"Join","text":"join_product(base, r) constructs a ProductManifold\n\nConstruct rank-r join parameter manifold as a plain ProductManifold. For Manifolds.Segre, uses flattened (Euclidean(1), Sphere, ..., Sphere) factors per component to preserve current CP parameter layout.\n\n\n\n\n\n","title":"TensorKitchen.join_product"},{"category":"function","location":"join/#TensorKitchen.SegreProduct","page":"Join","text":"SegreProduct(dims, r)\n\nProduct manifold Manifolds.Segre(dims) × ... × Manifolds.Segre(dims) (r factors). Each component point uses the Manifolds.Segre layout [[λ], x₁, …, x_d].\n\n\n\n\n\n","title":"TensorKitchen.SegreProduct"},{"category":"section","location":"btd/#Block-Term-Decomposition","page":"BTD","text":"A block term decomposition (BTD) with r blocks writes\n\nhat A = sum_i=1^r A_i\n\nwhere each block A_i is represented as a Tucker decomposition. At present, only homogeneous BTDs are supported, that is, all blocks must have the same multilinear rank.\n\nTo compute a block term decomposition of A with 10 blocks, each of multilinear rank (5, 4, 3), use\n\njulia> r = 10\njulia> mlrank = (5, 4, 3)\njulia> btd_res = btd(A, r, mlrank)\nBTDResult{Float64}\n Blocks: 10\n Rel. error: 0.2551559591470521\n\nThe blocks of btd_res can be obtained as follows:\n\nblocks = blocks(btd_res)\n\nEach block is represented as a Tucker decomposition, so we can access its core and factor matrices via:\n\nblk = blocks[1]\ncore(blk)\nfactors(blk)","title":"Block Term Decomposition"},{"category":"section","location":"btd/#BTD-Docs","page":"BTD","text":"","title":"BTD Docs"},{"category":"function","location":"btd/#TensorKitchen.btd","page":"BTD","text":"btd(A, blocks, ranks; kwargs...) returns a BTDResult\n\nComputes a block-term decomposition of A with blocks Tucker blocks, each with multilinear rank ranks. The solver first finds an initial point, then refines it. Returns a BTDResult.\n\nMain Options\n\ninit = :auto: Sets the algorithm to find the initial point. Possible options are:\n:auto: Uses a default BTD initializer. For solver = :als, this uses BTDHOSVDMultistartInit; otherwise, it uses an ALS warm start.\n:alswarm: Runs ALS first and uses the result as the initial point for manifold solver refinement.\ncustom initializer objects, e.g. BTDHOSVDMultistartInit(...).\nsolver = :rgd: Sets the algorithm for refinement. Possible options are:\n:rgd (default): Riemannian gradient descent.\n:als: Alternating least squares.\n:rcg: Riemannian conjugate gradient.\n:lbfgs: Limited-memory quasi-Newton refinement.\n\nExtended Options\n\ninit_point = nothing: Explicit initial point. If provided, it overrides the default initial point.\nwarm_init = BTDHOSVDMultistartInit(...): Searches for initial points using HOSVD for :alswarm, optionally screens them with short ALS runs, and returns the lowest-cost candidate.\nwarm_steps = 200: Once finding the best initial point, it runs this many ALS iterations to refine the initial point.\nwarm_block_method = :hooi or :sthosvd: Block update method used during warm start.\nwarm_block_maxiter = 20: Maximum number of inner iterations for each block update during warm start.\nwarm_rel_error_gate = 5e-2: Skips manifold refinement if the warm-start error is above this threshold.\nmaxiter = 500: Maximum number of Riemannian gradient descent iterations.\nstepsize = 0.01: Initial step size for line search in Riemannian gradient descent.\ntol = 1e-6: Convergence tolerance.\ngradient_mode = :riemannian: rgrad can be directly applied for manifold solvers. \nIf the model has a direct rgrad, it uses that.\nOtherwise it computes egrad and projects it to the tangent space.\nThis behavior is in src/solvers/abstract.jl (line 289).\nverbose = true: Enables progress output.\nblock_method = :hooi or :sthosvd: Block update method used for manifold solvers.\nblock_maxiter = 30: Maximum number of inner block-update iterations for manifold solvers.\nbtd_als_polish_maxiter = nothing: Number of final ALS polishing iterations. If nothing, an automatic budget is selected.\nThese settings are a robustness/quality feature for BTD-ALS. They are not used for manifold solvers.\nmax_stagnation_restarts = 1: More retries after a bad stagnated ALS pass. Higher values can improve solution quality, but increases runtime.\nstagnation_rel_error = 1e-4: This is a “bad final error” cutoff, not an improvement threshold. Lower value means restarts trigger more easily.\nrestart_candidates = 24: It usually improves chance of finding a better basin, but cost grows roughly linearly.\nrestart_screening_steps = 5: It uses this many quick ALS steps to screen restart candidates.\nrestart_block_maxiter = 20: Inner block-update limit during restart screening for BTD-ALS.\nrestart_seed = nothing: Optional random seed for restart generation for BTD-ALS.\n\nExample\n\njulia> using Random\njulia> Random.seed!(0)\njulia> A = randn(20, 15, 10); blocks = 10; ranks = (5, 4, 3)\njulia> res = btd(A, blocks, ranks; verbose = false)\nBTDResult{Float64}\n Blocks: 10\n Rel. error: 0.2625821087015455\n\n\n\n\n\n","title":"TensorKitchen.btd"},{"category":"type","location":"btd/#TensorKitchen.BTDResult","page":"BTD","text":"BTDResult{T}\n\nResult of block-term decomposition (btd); block components expose Tucker structure through accessors like core(blk), factors(blk), and blk.tensor.\n\nsolver_info: solver-specific diagnostics/metadata (NamedTuple). Typical keys include BTD-ALS restart diagnostics (total_iterations, stagnation_restarts, restart_rel_error_history) and BTD-TSD run settings (schedule, block_repeats, block_count, stepsize).\n\n\n\n\n\n","title":"TensorKitchen.BTDResult"},{"category":"function","location":"btd/#TensorKitchen.blocks","page":"BTD","text":"blocks(r::BTDResult)\n\nReturn the Tucker block components of a block-term decomposition result.\n\n\n\n\n\n","title":"TensorKitchen.blocks"},{"category":"method","location":"btd/#TensorKitchen.reconstruct-Tuple{BTDResult}","page":"BTD","text":"reconstruct(td::TuckerResult) reconstructs the tensor from Tucker decomposition\n\nA = S ×₁ U₁ ×₂ U₂ ⋯ ×_d U_d\n\n\n\n\n\n","title":"TensorKitchen.reconstruct"},{"category":"section","location":"tucker/#Tucker-Decomposition","page":"Tucker","text":"Approximating A by a Tucker decomposition\n\nhat A = C times_1 U times_2 V times_3 W\n\nwith multilinear rank mlrank can be computed as follows.\n\njulia> mlrank = (5, 4, 3)\njulia> tucker_res = tucker(A, mlrank)\nTuckerResult{Float64, 3}\n Original size: (20, 15, 10)\n Core size: (5, 4, 3)\n Multilinear rank: (5, 4, 3)\n Compression: 12.0x\n\nThe core C and the factor matrices (U V W) of the decomposition can be accessed as follows.\n\ncore(tucker_res)\nfactors(tucker_res)","title":"Tucker Decomposition"},{"category":"section","location":"tucker/#Tucker-Docs","page":"Tucker","text":"","title":"Tucker Docs"},{"category":"function","location":"tucker/#TensorKitchen.tucker","page":"Tucker","text":"tucker(A, ranks; method = :sthosvd, kwargs...) returns a TuckerResult\n\nComputes a Tucker decomposition of A with multilinear rank ranks.\n\nMain Options\n\nmethod = :sthosvd: Sets the Tucker decomposition algorithm. Possible options are:\n:sthosvd (default): Sequentially Truncated HOSVD. A direct one-pass decomposition, mainly used as a fast standalone Tucker approximation or as the default initializer for :hooiFast, deterministic, and usually a good initial point.\n:hooi: High-Order Orthogonal Iteration. Iteratively refines the Tucker factors, initialized by ST-HOSVD by default.\n\nExtended Options\n\nFor method = :hooi:\n\nmaxiter = 50: Maximum number of HOOI iterations.\ntol = 1e-8: Convergence tolerance based on change in relative reconstruction error.\ninit = :sthosvd\n:sthosvd: Uses ST-HOSVD to initialize the Tucker factors.\nTuckerResult: Uses an existing Tucker decomposition as the initial point.\n\nExample\n\njulia> using Random\njulia> Random.seed!(0)\njulia> A = randn(20, 15, 10); ranks = (5, 4, 3)\njulia> res = tucker(A, ranks; verbose = false)\nTuckerResult{Float64, 3}\n Original size: (20, 15, 10)\n Core size: (5, 4, 3)\n Multilinear rank: (5, 4, 3)\n Compression: 12.0x\n\n\nThe core tensor and factor matrices can be accessed by\n\ncore(res)\nfactors(res)\n\nA tensor approximation can be reconstructed by\n\nreconstruct(res)\n\nor equivalently\n\nreconstruct_tucker(core(res), factors(res))\n\nNotes\n\n:sthosvd is not an iterative solver. It directly returns a TuckerResult and does not expose solver-style outputs such as iteration counts or convergence diagnostics.\n:hooi is the iterative refinement method in the current Tucker implementation.\nImportant distinction: For the current Tucker implementation, do not use solver = :rgd, init = :auto, or manifold solvers like RGD, RCG, LBFGS, etc.\n\n\n\n\n\n","title":"TensorKitchen.tucker"},{"category":"type","location":"tucker/#TensorKitchen.TuckerResult","page":"Tucker","text":"TuckerResult{T, N}\n\nStores a Tucker decomposition: core tensor and factor matrices.\n\ncore::Array{T,N} — core tensor\nfactors::Vector{Matrix{T}} — orthonormal factor matrices\nprocessing_order::Vector{Int} — order modes were processed\nsingular_values::Vector{Vector{T}} — singular values per truncation\n\n\n\n\n\n","title":"TensorKitchen.TuckerResult"},{"category":"method","location":"tucker/#TensorKitchen.core-Tuple{TuckerResult}","page":"Tucker","text":"core(td::TuckerResult)\n\nReturn the Tucker core tensor for a Tucker result.\n\n\n\n\n\n","title":"TensorKitchen.core"},{"category":"method","location":"tucker/#TensorKitchen.factors-Tuple{TuckerResult}","page":"Tucker","text":"factors(td::TuckerResult)\n\nReturn the Tucker factor matrices.\n\n\n\n\n\n","title":"TensorKitchen.factors"},{"category":"method","location":"tucker/#TensorKitchen.multilinear_rank-Tuple{TuckerResult}","page":"Tucker","text":"multilinear_rank(td::TuckerResult)\n\nReturn the Tucker multilinear rank tuple, i.e. the size of the core tensor.\n\n\n\n\n\n","title":"TensorKitchen.multilinear_rank"},{"category":"method","location":"tucker/#TensorKitchen.factor_dims-Tuple{TuckerResult}","page":"Tucker","text":"factor_dims(td::TuckerResult)\n\nReturn the original mode dimensions represented by the Tucker factor matrices.\n\n\n\n\n\n","title":"TensorKitchen.factor_dims"},{"category":"method","location":"tucker/#TensorKitchen.reconstruct-Tuple{TuckerResult}","page":"Tucker","text":"reconstruct(td::TuckerResult) reconstructs the tensor from Tucker decomposition\n\nA = S ×₁ U₁ ×₂ U₂ ⋯ ×_d U_d\n\n\n\n\n\n","title":"TensorKitchen.reconstruct"},{"category":"section","location":"cpd/#CPD","page":"CPD","text":"Here is how to approximate a tensor A by a CPD of rank r.\n\njulia> using TensorKitchen\njulia> A = randn(20, 15, 10)\njulia> r = 35\njulia> res = cpd(A, r)\nCPDResult{Float64}\n Order: 3\n Dimensions: (20, 15, 10)\n Rank: 35\n Rel. error: 0.4359141301703327\n\nNow, res contains a CP approximation of the 3-way tensor A,\n\nhat A = sum_i=1^r lambda_i a_i otimes b_i otimes c_i\n\nIt approximates A with relative error about 0.436.\n\nWe access the decomposition as follows.\n\nλ = weights(res)\nU = factors(res)\n\nHere, U is a triple of matrices (ABC), where the columns of A are the a_i and so on. These are called factor matrices.\n\nWe get the whole reconstructed tensor by \n\n = reconstruct(res)","title":"CPD"},{"category":"section","location":"cpd/#CPD-Docs","page":"CPD","text":"","title":"CPD Docs"},{"category":"function","location":"cpd/#TensorKitchen.cpd","page":"CPD","text":"cpd(A, r; kwargs...)\n\nComputes a rank-r CP approximation of A in two steps: (1) the first step finds an initial point; (2) the second step refines the initial point. Returns a CPDResult. If r is omitted, uses the smallest tensor mode as a heuristic rank.\n\nMain Options\n\ninit = :auto: Sets the algorithm to find the initial point. Possible options are:\n:auto: Uses a default CPD initializer. For solver = :als, this uses TuckerInit; otherwise, it uses an ALS warm start.\n:alswarm: Runs ALS first and uses the result as the initial point for refinement.\ncustomized initial point:\n:tucker (default when solver = :als): Uses a default Tucker initializer.\n:random: Uses a random initial point.\n:hosvd: Uses a HOSVD initial point.\nsolver = :rgd: Sets the algorithm for refinement. Possible options are:\nrgd (default): Riemannian gradient descent\nrgd_fixed: Riemannian gradient descent with fixed step size\nrcg: Riemannian conjugate gradient\nals: Alternating Least Squares\n\nExtended Options\n\np0 = nothing: Explicit initial point. If provided, it overrides the default initial point.\n:alswarm: ALS warm start option.\nwarm_init = TuckerInit(): Before finding the warm start initial point, this sets the good starting point for ALS.\nwarm_steps = 500: Once finding the best initial point from warm_init, it runs this many ALS iterations to refine the initial point.\nmaxiter = 500: Maximum number of Riemannian gradient descent iterations.\nstepsize = 1.0: Initial step size for line search in Riemannian gradient descent.\ntol = 1e-6: Convergence tolerance.\ngradient_mode = :riemannian: Gradient rule for manifold solvers. \nIf the model has a direct rgrad, it uses that.\nOtherwise it computes egrad and projects it to the tangent space.\nThis behavior is in src/solvers/abstract.jl (line 289).\ngeometry = :canonical: Sets the geometry of the manifold. Possible options are:\n:canonical: Standard CPD parameterization with the usual Euclidean factors and canonical Riemannian gradient handling. Best default for general unconstrained CPD.\n:squaring_metric: Nonnegative geometry based on squared latent coordinates. Enforces nonnegativity indirectly, but can become ill-conditioned near zero.\n:softplus_metric: Nonnegative geometry uses a regularized pullback-inspired geometry induced by the softplus chart. Smoother and usually more stable near zero than :squaring_metric.\n:native: Native CP manifold geometry using the model’s intrinsic CP/Segre representation not for nonnegative=true. Best for structured join layouts with Manifolds.Segre summands.\nverbose = true: Enables progress output.\nnonnegative::Bool = false: Nonnegative CPD option to be selected by the user. (same as nncpd)\npullback_eps = 1e-8: Regularization parameter for pullback-style nonnegative geometries.\n\nNotes\n\nsolver = :als does not use manifold geometry. In that case:\ngeometry must be :canonical\ngradient_mode is ignored except for validation\n:squaring_metric and :softplus_metric require nonnegative = true.\nWhen nonnegative = true, cpd(...) routes to nncpd(...). In that route:\nif solver != :als and geometry is left at :canonical, the effective geometry becomes :softplus_metric\nif stepsize is left at 1.0, the effective default becomes 0.01\nif init = :tucker, the effective initializer becomes :alswarm\n\nExample\n\njulia> A = randn(20, 15, 10); r = 35\njulia> res = cpd(A, r)\nCPDResult{Float64}\n Order: 3\n Dimensions: (20, 15, 10)\n Rank: 35\n Rel. error: 0.4359141301703327\n\n\n\n\n\n","title":"TensorKitchen.cpd"},{"category":"function","location":"cpd/#TensorKitchen.nncpd","page":"CPD","text":" nncpd(A, r; kwargs...)\n\nComputes a nonnegative rank-r CP approximation of A in two steps: (1) the first step finds an initial point; (2) the second step refines the initial point. Returns a CPDResult. If r is omitted, uses the smallest tensor mode as a heuristic rank. cpd(A, r; nonnegative=true, ...) routes here and adopts the same effective defaults.\n\nOptions\n\nThe options are the same as for cpd.\n\nGeometry guide:\n\ngeometry=:softplus_metric Default and usually the safest choice. \ngeometry=:squaring_metric Uses a regularized pullback-inspired geometry induced by the squaring chart.\ngeometry=:canonical Plain nonnegative CP coordinates without the pullback-style manifold geometry. This is the natural choice with solver=:als.\n\nExample\n\njulia> A = randn(20, 15, 10); r = 35;\njulia> B = abs.(A)\njulia> nncpd(B, r)\nCPDResult{Float64}\n Order: 3\n Dimensions: (20, 15, 10)\n Rank: 35\n Rel. error: 0.3765605093526155\n\n\n\n\n\n","title":"TensorKitchen.nncpd"},{"category":"type","location":"cpd/#TensorKitchen.CPDResult","page":"CPD","text":"CPDResult{T}\n\nResult of a Canonical Polyadic Decomposition.\n\nStores the decoded CP representation together with solver diagnostics:\n\ncomponents: rank-one tensor components\nweights: component weights\nfactors: factor matrices\ncost: final objective function value at the returned solution\nrel_error: final relative reconstruction error\ngrad_norm: norm of the final optimization gradient reported by the solver; for manifold solvers this is the Riemannian gradient norm\niterations: number of refinement iterations\nconverged: whether the solver reported convergence\nsolver: solver optimization method used to produce the result\nsolver_info: solver-specific diagnostics/metadata (NamedTuple). Typical keys include:\ninitial_stepsize_eff (RGD), memory_size (LBFGS), cautious_update (LBFGS),\ninitial_scale, linesearch, has_preconditioner (LBFGS), and nncp_pullback_eps (NNCP).\n\n\n\n\n\n","title":"TensorKitchen.CPDResult"},{"category":"method","location":"cpd/#TensorKitchen.weights-Tuple{CPDResult}","page":"CPD","text":"weights(r::CPDResult)\n\nReturn the CP component weights stored in a CPD result.\n\n\n\n\n\n","title":"TensorKitchen.weights"},{"category":"method","location":"cpd/#TensorKitchen.factors-Tuple{CPDResult}","page":"CPD","text":"factors(res::CPDResult)\n\nReturn the CP factor matrices of res as a vector [U₁, U₂, ..., U_N], where each U_m has size size(A, m) × rank.\n\n\n\n\n\n","title":"TensorKitchen.factors"},{"category":"method","location":"cpd/#TensorKitchen.reconstruct-Tuple{CPDResult}","page":"CPD","text":"reconstruct(td::TuckerResult) reconstructs the tensor from Tucker decomposition\n\nA = S ×₁ U₁ ×₂ U₂ ⋯ ×_d U_d\n\n\n\n\n\n","title":"TensorKitchen.reconstruct"},{"category":"section","location":"#TensorKitchen.jl-Documentation","page":"Home","text":"TensorKitchen.jl is a Julia package for tensor decompositions.\n\nCPD documentation\nTucker documentation\nBTD documentation\nJoin decomposition documentation\nUtilities\nPipeline\nReferences","title":"TensorKitchen.jl Documentation"},{"category":"section","location":"#Notes","page":"Home","text":"The package is currently at a pre-alpha stage. \n\nThe implementation is based on combining algebraic algorithms like ALS (see, e.g., the textbook by Kolda and Ballard) and Riemannian optimization from Manopt.jl.\n\nWhat currently works is \n\nCanonical Polyadic Decomposition (CPD)\nTucker Decomposition\nNonnegative Canonical Polyadic Decomposition (NNCPD)\nBlock Term Decomposition (BTD)\nJoin Decompositions\n\n\n\nThe next updates will include \n\nHandling of swamps/plateaus in the optimization step\nDocumentation\nImproved User Interface\nGPU Support \nLL1 Decomposition (3-way specialized BTD)\nSymmetric CP / Waring Decomposition\nPartially Symmetric CP\nTensor Trains","title":"Notes"},{"category":"section","location":"utils/#General-Utilities","page":"Utilities","text":"","title":"General Utilities"},{"category":"function","location":"utils/#TensorKitchen.save_result","page":"Utilities","text":"save_result(path::AbstractString, result)\n\nSave a result to a file.\n\n\n\n\n\n","title":"TensorKitchen.save_result"},{"category":"function","location":"utils/#TensorKitchen.load_result","page":"Utilities","text":"load_result(path::AbstractString)\n\nLoad a result from a file\n\n\n\n\n\n","title":"TensorKitchen.load_result"}] +[{"category":"section","location":"references/#References","page":"References","text":"","title":"References"},{"category":"section","location":"references/#General-Tensor-Decomposition","page":"References","text":"Tensor decompositions (CP, Tucker): T. G. Kolda and B. W. Bader, \"Tensor decompositions and applications,\" SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.","title":"General Tensor Decomposition"},{"category":"section","location":"references/#Tucker-Methods","page":"References","text":"HOSVD: L. De Lathauwer, B. De Moor, and J. Vandewalle, \"A multilinear singular value decomposition,\" SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1253–1278, 2000.\nST-HOSVD: N. Vannieuwenhoven, R. Vandebril, K. Meerbergen, \"A new truncation strategy for the higher-order singular value decomposition,\" SIAM J. Sci. Comput., vol. 34, no. 2, pp. A1027–A1052, 2012.\nHOOI: L. De Lathauwer, B. De Moor, and J. Vandewalle, \"On the best rank-1 and rank-(R1,R2,...,R_N) approximation of higher-order tensors,\" SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1324–1342, 2000.","title":"Tucker Methods"},{"category":"section","location":"references/#Block-and-Structured-Models-(BTD-/-LL1)","page":"References","text":"Block-term decomposition (BTD): L. De Lathauwer, \"Decompositions of a higher-order tensor in block terms—Part I: Lemmas for partitioned matrices,\" SIAM J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1022–1032, 2008.\nL. De Lathauwer, \"Decompositions of a higher-order tensor in block terms—Part II: Definitions and uniqueness,\" SIAM J. Matrix Anal. Appl., vol. 30, no. 3, pp. 1033–1066, 2008.\nBTD-ALS: L. De Lathauwer and D. Nion, \"Decompositions of a higher-order tensor in block terms—Part III: Alternating least squares algorithms,\" SIAM Journal on Matrix Analysis and Applications, vol. 30, no. 3, pp. 1067–1083, 2008. PDF.","title":"Block and Structured Models (BTD / LL1)"},{"category":"section","location":"references/#Join-decompositions","page":"References","text":"Conditioning of join decompositions: P. Breiding and N. Vannieuwenhoven, \"The condition number of join decompositions,\" SIAM Journal on Matrix Analysis and Applications, vol. 39, no. 1, pp. 287–309, 2018. arXiv:1611.08117 (PDF).","title":"Join decompositions"},{"category":"section","location":"references/#Riemannian-Optimization-and-Julia-Ecosystem","page":"References","text":"Riemannian trust-region / Gauss–Newton for canonical rank (CP) approximation: P. Breiding and N. Vannieuwenhoven, \"A Riemannian Trust Region Method for the Canonical Tensor Rank Approximation Problem,\" SIAM Journal on Optimization, vol. 28, no. 3, pp. 2435–2465, 2018. arXiv:1709.00033 (PDF).\nRiemannian optimization: P.-A. Absil, R. Mahony, and R. Sepulchre, Optimization Algorithms on Matrix Manifolds. Princeton University Press, 2008.\nJulia manifold optimization ecosystem: R. Bergmann et al., ManifoldsBase.jl, Manifolds.jl, and Manopt.jl.","title":"Riemannian Optimization and Julia Ecosystem"},{"category":"section","location":"PIPELINE/#TensorKitchen-Pipeline","page":"Pipeline","text":"This document explains how public APIs route into models, solvers, and result converters.","title":"TensorKitchen Pipeline"},{"category":"section","location":"PIPELINE/#Public-entry-points","page":"Pipeline","text":"CP Decomposition cpd(A, r; ...) \nNonnegative CP Decomposition nncpd(A, r; ...) \nBlock Term Decomposition btd(A, blocks, ranks; ...) \nTucker Decomposition tucker(A, ranks; method=...) \nJoin Decomposition approx(...) ","title":"Public entry points"},{"category":"section","location":"PIPELINE/#Default-behavior-(quick-reference)","page":"Pipeline","text":"cpd(A, r):\ninit = :alswarm\nsolver = :rgd\nnncpd(A, r):\ninit = :alswarm\nsolver = :rgd\nbtd(A, blocks, ranks):\ninit = :alswarm\nwarm_steps = 200\nwarm_init = BTDHOSVDMultistartInit(candidates=64, screening_steps=10, block_maxiter=12)\nwarm_rel_error_gate = nothing (run manifold refinement by default; set e.g. 5e-2 to short-circuit on poor warm starts)\nsolver = :rgd\nfinal BTD-ALS polish enabled by default for non-ALS solvers\nmax_stagnation_restarts = 1 (retry with stronger multistart when ALS fit-change stalls at high rel-error)\ntucker(A, ranks):\nmethod = :sthosvd\napprox(model::JoinModel):\ninit = :random\nsolver = :rgd","title":"Default behavior (quick reference)"},{"category":"section","location":"PIPELINE/#Core-execution-architecture","page":"Pipeline","text":"Most optimization APIs share this core pattern:\n\nBuild a model (JoinModel + backend)\nCall _solve_model(...)\nConvert to a public result struct\n\n_solve_model lives in src/solvers/solve_dispatch.jl and is the common symbol-to-solver dispatch layer (:rgd, :rcg, :lbfgs, :als, :btd_tsd).","title":"Core execution architecture"},{"category":"section","location":"PIPELINE/#API-flows","page":"Pipeline","text":"","title":"API flows"},{"category":"section","location":"PIPELINE/#CPD-(cpd,-nncpd)","page":"Pipeline","text":"cpd(A, r; ...):\n\nBuild JoinModel(A, r; geometry=...) with CPDBackend\nNormalize/validate options (solver, geometry, gradient_mode, normalization policy)\nSolve through _solve_model(...)\nConvert to CPDResult\n\nNotes:\n\n:als means CP-ALS.\nManifold solvers (:rgd, :rgd_fixed, :rcg, :lbfgs) share dispatch with other pipelines.\nFor solver != :als, init = :auto resolves to :alswarm, so CPD and NNCPD start from an ALS warm point before manifold refinement.\nGeneric approx(...) does not use CPD's ALS warm-start path unless it auto-routes to cpd(...).","title":"CPD (cpd, nncpd)"},{"category":"section","location":"PIPELINE/#BTD-(btd)","page":"Pipeline","text":"btd(A, blocks, ranks; ...):\n\nBuild a uniform Tucker family via TuckerJoin(...)\nWrap as JoinModel with BTDBackend\nChoose effective initializer:\nsolver == :als: use requested init directly (default multistart)\nsolver != :als: use BTDALSWarmStartInit(...) so first-order methods start from a good BTD-ALS warm point\nIf the warm-start rel-error exceeds warm_rel_error_gate, return the warm BTD-ALS result directly\nOtherwise solve through _solve_model(...)\nIf solver != :als, optionally polish with BTD-ALS (btd_als_polish_maxiter)\nConvert to BTDResult\n\nPolish step usefulness:\n\nUsually helpful for a small final rel_error reduction after RGD converges near a good basin.\nMost useful for quality-focused runs (benchmarks, final fits).\nCan be skipped for speed-sensitive runs (btd_als_polish_maxiter=0) when small extra gains are not worth runtime.\n\nBTD-specific initialization options:\n\n:hosvd: sequential block initialization on residual\n:hosvd_multistart: HOSVD subspace split candidates, optional screening ALS, keep lowest-cost candidate\n:alswarm: short BTD-ALS warm-start wrapper around base initializer\n\nBTD-ALS stabilization behavior:\n\nTracks per-iteration fit change (|rel_t - rel_{t-1}|)\nDetects stagnation when fit change is tiny but rel_error remains high\nCan restart from fresh multistart pool (max_stagnation_restarts)\nReports true final Riemannian gradient norm (grad_norm) ","title":"BTD (btd)"},{"category":"section","location":"PIPELINE/#Tucker-(tucker)","page":"Pipeline","text":"tucker(A, ranks; method=...) does not use _solve_model. It dispatches directly to decomposition routines:\n\n:sthosvd\n:hooi","title":"Tucker (tucker)"},{"category":"section","location":"PIPELINE/#Generic-approx(...)-routing","page":"Pipeline","text":"approx(manifolds, target; dispatch=:auto) routes by manifold family:\n\nuniform Manifolds.Segre -> cpd(...)\nuniform Manifolds.Tucker matching target shape/rank -> btd(...)\nmixed or non-uniform family -> generic JoinModel(...) path -> ApproxResult\n\ndispatch=:cpd, :btd, and :generic force behavior.\n\nFor the generic JoinModel path, approx(...) starts from init = :random by default and then runs the selected manifold solver. It does not run an ALS warm-start stage, because a general join component does not necessarily expose factor matrices or least-squares block updates.","title":"Generic approx(...) routing"},{"category":"section","location":"PIPELINE/#Result-types-and-post-processing","page":"Pipeline","text":"CPDResult\nBTDResult\nTuckerResult\nApproxResult\n\nCommon utilities:\n\nreconstruct(result)\nrel_error(A, result)","title":"Result types and post-processing"},{"category":"section","location":"PIPELINE/#File-map","page":"Pipeline","text":"API entry points: src/api/approx.jl, src/api/cpd.jl, src/api/nncpd.jl, src/api/btd.jl, src/api/tucker.jl","title":"File map"},{"category":"section","location":"join/#Join-Decomposition","page":"Join","text":"A Join Decomposition of a vector xinmathbb R^N is a decomposition of the form x = x_1+cdots+x_r, where x_iin M_i and M_isubset mathbb R^N is a given embedded manifold. \n\nFor instance, we can approximate a point p = (12 04)inmathbb R^2 by x=x_1+x_2, where x_1x_2in S^1 are points on the circle:\n\njulia> using Manifolds\njulia> p = [1.2, 0.4]\njulia> S = Sphere(1)\njulia> join_res = approx((S, S), p)\nApproxResult{Float64}\n Components: 2\n Rel. error: 0.0007114699550529457\n\nWe access the decomposition as follows.\n\ncomponents(join_res)\nreconstruct(join_res)\n\n
","title":"Join Decomposition"},{"category":"section","location":"join/#Generic-Join-Approximation","page":"Join","text":"approx(...) is the main frontend for join decomposition. It works in two stages:\n\nbuild an initial point\nrefine it with the selected solver","title":"Generic Join Approximation"},{"category":"section","location":"join/#Supported-Forms","page":"Join","text":"approx(model; kwargs...) Solve an already constructed JoinModel.\napprox(manifolds, target; kwargs...) Build a join from a tuple/vector of component manifolds.\napprox(M::ProductManifold, target; kwargs...) Use the product-manifold factors as join components.\napprox(base, r, target; kwargs...) Repeat a single base manifold r times to build a join.\napprox(base, target; kwargs...) Build a single-component join.\n\nExamples:\n\njulia> target = [1.2, 0.4, -0.3]\njulia> model = JoinModel(Manifolds.Sphere(2), target)\njulia> approx(model; maxiter = 50, verbose = false)\nApproxResult{Float64}\n Components: 1\n Rel. error: 0.23076923076923075\n\njulia> target = randn(2, 3)\njulia> approx((Manifolds.Segre((2, 3)), Manifolds.Segre((2, 3))), target; verbose = false)\nCPDResult{Float64}\n\njulia> target = [1.2, 0.4, -0.3]\njulia> approx(Manifolds.Sphere(2), 2, target; verbose = false)\nApproxResult{Float64}","title":"Supported Forms"},{"category":"section","location":"join/#Routing","page":"Join","text":"With dispatch = :auto:\n\nuniform Manifolds.Segre components route to cpd(...)\nuniform Manifolds.Tucker components route to btd(...)\notherwise the generic JoinModel(...) path is used and ApproxResult is returned\n\nYou can also force the route explicitly:\n\ndispatch = :generic\ndispatch = :cpd\ndispatch = :btd\n\nForced routing is validated:\n\ndispatch = :cpd requires all components to be Manifolds.Segre with identical factor_dims\ndispatch = :btd requires all components to be Manifolds.Tucker with identical factor_dims and compatible multilinear rank","title":"Routing"},{"category":"section","location":"join/#Generic-Join-Options","page":"Join","text":"For the generic join path:\n\ninit = :random Default initializer. Built-in initializer support depends on the component manifolds.\nsolver = :rgd Supported generic-join solver options are:\n:rgd\n:rgd_fixed\n:rcg\n:lbfgs\n\nOther common options:\n\np0 = nothing\nmaxiter = 500\nstepsize = 1.0\ntol = 1e-6\ngradient_mode = :riemannian\nverbose = true\nvector_transport_method = nothing\n\nBuilt-in init support by component family:\n\nSphere: :random, :deterministic, :target\nSegre: :random, :deterministic\nTucker: :random, :tucker, :tucker_diag, :sthosvd\nother manifolds: :random","title":"Generic Join Options"},{"category":"section","location":"join/#Notes","page":"Join","text":"Generic joins require every component manifold to embed into the same flattened ambient length as target.\nsolver = :als is not available for a truly generic JoinModel. ALS may still be available when approx(...) auto-routes to cpd(...) or btd(...).","title":"Notes"},{"category":"section","location":"join/#Join-Decomposition-Docs","page":"Join","text":"","title":"Join Decomposition Docs"},{"category":"function","location":"join/#TensorKitchen.approx","page":"Join","text":"Generic Join Approximation\n\napprox(...) is the main frontend for join decomposition, which works in two stages:\n\nbuild an initial point\nrefine it with the selected solver\n\nSupported Forms\n\napprox(model; kwargs...) : It is for an already constructed join model (JoinModel(...)) and routes to the generic join solver. Model can be a JoinModel of a tuple of manifolds, a ProductManifold, or a single manifold.\nmodel means an already constructed JoinModel.\nIt fully fixes the decomposition structure and target.\napprox(model; ...) just solves that model.\nExample: If you want to approximate a point on the sphere, you can build a JoinModel and then use approx to refine it.\n\ntarget = [1.2, 0.4, -0.3]\nmodel = JoinModel(Manifolds.Sphere(2), target)\napprox(model; maxiter = 100, verbose = false)\n# returns an ApproxResult\n\napprox(manifolds, target; kwargs...) : Builds a generic join model from a tuple of existing manifolds and routes to according to dispatch. Example:\n\ntarget = randn(2, 3)\napprox((Manifolds.Segre((2, 3)), Manifolds.Segre((2, 3))), target; verbose = false)\n# This builds a join with 2 copies of Manifolds.Segre((2, 3))\".\n\napprox(M::ProductManifold, target; kwargs...) : Uses the factors of a product manifold and route to CPD, BTD, or the generic join solver according to dispatch. \nM::ProductManifold means that you already have a product manifold whose factors are the join components.\napprox(M, target; ...) uses those factors directly.\nIt is more explicit than base, because the components are already listed.\nExample:\n\n# Example 1\ntarget = randn(2, 3)\nM = ProductManifold(Manifolds.Segre((2, 3)), Manifolds.Segre((2, 3)))\napprox(M, target; verbose = false)\n# returns an CPDResult\n\n# Example 2\ntarget = randn(4, 3, 2)\nM = ProductManifold(\n Manifolds.Tucker((4, 3, 2), (2, 2, 2)),\n Manifolds.Tucker((4, 3, 2), (2, 2, 2)),\n)\napprox(M, target; verbose = false)\n# returns an BTDResult\n\napprox(base, r, target; kwargs...) : builds a rank-r Segre join and routes to CPD when base isa Manifolds.Segre and BTD when base isa Manifolds.Tucker unless generic\n\ndispatch is explicitly requested. - base means one manifold template, not yet a full join. - approx(base, r, target; ...) repeats that same manifold r times to build a join. - approx(base, target; ...) builds a one-component join.\n\ntarget = randn(2, 3)\napprox(Manifolds.Segre((2, 3)), 2, target; verbose = false)\n# returns an CPDResult\n\napprox(base, target; kwargs...) : Builds a single-component generic join and route to the generic join solver. Example:\n\ntarget = [1.2, 0.4, -0.3] \napprox(Manifolds.Sphere(2), target; verbose = false)\n# returns an ApproxResult\n\nBy default, approx auto-routes by manifold family:\nuniform Manifolds.Segre summands calls cpd(...)\nuniform Manifolds.Tucker summands calls btd(...)\notherwise calls JoinModel(...) and returns a ApproxResult\n\nReturn Types\n\nDepending on the manifold family, approx(...) may return:\nApproxResult for the generic join path\nCPDResult when auto-routed to cpd(...)\nBTDResult when auto-routed to btd(...)\n\nMain Options\n\nFor the generic join path:\n\ninit = :random: Sets the algorithm to find the initial point.\nsolver = :rgd: Sets the algorithm for refinement. Possible options are:\nrgd (default): Riemannian gradient descent\nrgd_fixed: Riemannian gradient descent with fixed step size\nrcg: Riemannian conjugate gradient\nlbfgs: Limited-memory quasi-Newton\n\n##Notes##\n\n:als is not a solver option for approx(...). However, if approx(...) auto-routes to cpd(...) or btd(...), then those specialized pipelines may support ALS separately.\nwarm_steps and warm_init are not part of the generic approx(...) path. Generic joins start from random initial point and then use manifold solvers for refinement.\nFor generic mixed joins, use manifold solvers such as :rgd, :rcg, or :lbfgs.\n\n\n\n\n\napprox(M::ProductManifold, target; dispatch=:auto, kwargs...)\n\nUse the factors of a product manifold as join components and route to CPD, BTD, or the generic join solver according to dispatch.\n\n\n\n\n\napprox(base::Manifolds.Segre, r, target; dispatch=:auto, kwargs...)\n\nBuild a rank-r Segre join and route to the CPD pipeline unless generic dispatch is explicitly requested.\n\n\n\n\n\napprox(base::Manifolds.Tucker, r, target; dispatch=:auto, kwargs...)\n\nBuild a r-block Tucker join and route to the BTD pipeline unless generic dispatch is explicitly requested.\n\n\n\n\n\napprox(base::AbstractManifold, r, target; dispatch=:auto, kwargs...)\n\nFallback rank-r join constructor for non-specialized manifolds. Forced CPD or BTD dispatch is rejected because the base manifold family is not known.\n\n\n\n\n\napprox(base::AbstractManifold, target; dispatch=:auto, kwargs...)\n\nSingle-component generic approximation fallback. Use this when no CPD/BTD family-specific routing is intended.\n\n\n\n\n\n","title":"TensorKitchen.approx"},{"category":"type","location":"join/#TensorKitchen.ApproxResult","page":"Join","text":"ApproxResult{T}\n\nGeneric result of approx(manifolds, target) (generic join decomposition).\n\npoint: final point on the join manifold\ncomponents: extracted component descriptions\ncost: the value of the optimization objective at the final returned point, that is typically the least-squares objective.\nrel_error: relative error of the decomposition which is the ratio of the cost to the target norm. is scale-normalized and easier to compare across problems\ngrad_norm: norm of the final optimization gradient reported by the solver; for manifold solvers this is typically the Riemannian gradient norm\niterations: number of iterations used\nconverged: whether the decomposition converged\nsolver: solver used for the decomposition\nsolver_info: solver-specific diagnostics/metadata (NamedTuple)\n\n\n\n\n\n","title":"TensorKitchen.ApproxResult"},{"category":"method","location":"join/#TensorKitchen.reconstruct-Tuple{ApproxResult}","page":"Join","text":"reconstruct(res::ApproxResult)\n\nReconstruct the dense ambient object represented by a generic join approximation result by summing its component tensors.\n\n\n\n\n\n","title":"TensorKitchen.reconstruct"},{"category":"function","location":"join/#TensorKitchen.join_product","page":"Join","text":"join_product(base, r)\n\nDecides how to expand base into r components:\n\nManifolds.Segre → flattened (Euclidean(1), Sphere, ...) × r (one λ + spheres per rank-1).\nManifolds.Tucker / ProductManifold → repeat each factor r times.\nGeneric manifold → ProductManifold(base, base, ..., base).\n\n\n\n\n\n","title":"TensorKitchen.join_product"},{"category":"function","location":"join/#TensorKitchen.SegreProduct","page":"Join","text":"SegreProduct(dims, r)\n\nProduct manifold Manifolds.Segre(dims) × ... × Manifolds.Segre(dims) (r factors).\nEach component point uses the Manifolds.Segre layout [[λ], x₁, …, x_d].\nSegreProduct is used to build a join model for CPD.\n\n\n\n\n\n","title":"TensorKitchen.SegreProduct"},{"category":"section","location":"btd/#Block-Term-Decomposition","page":"BTD","text":"A block term decomposition (BTD) with r blocks writes\n\nhat A = sum_i=1^r A_i\n\nwhere each block A_i is represented as a Tucker decomposition. At present, only homogeneous BTDs are supported, that is, all blocks must have the same multilinear rank.\n\nTo compute a block term decomposition of A with 10 blocks, each of multilinear rank (5, 4, 3), use\n\njulia> r = 10\njulia> mlrank = (5, 4, 3)\njulia> btd_res = btd(A, r, mlrank)\nBTDResult{Float64}\n Blocks: 10\n Rel. error: 0.2551559591470521\n\nThe blocks of btd_res can be obtained as follows:\n\nblocks = blocks(btd_res)\n\nEach block is represented as a Tucker decomposition, so we can access its core and factor matrices via:\n\nblk = blocks[1]\ncore(blk)\nfactors(blk)","title":"Block Term Decomposition"},{"category":"section","location":"btd/#BTD-Docs","page":"BTD","text":"","title":"BTD Docs"},{"category":"function","location":"btd/#TensorKitchen.btd","page":"BTD","text":"btd(A, blocks, ranks; kwargs...) returns a BTDResult\n\nComputes a block-term decomposition of A with blocks Tucker blocks, each with multilinear rank ranks. The solver first finds an initial point, then refines it. Returns a BTDResult.\n\nMain Options\n\ninit = :auto: Sets the algorithm to find the initial point. Possible options are:\n:auto: Uses a default BTD initializer. For solver = :als, this uses BTDHOSVDMultistartInit; otherwise, it uses an ALS warm start.\n:alswarm: Runs ALS first and uses the result as the initial point for manifold solver refinement.\ncustom initializer objects, e.g. BTDHOSVDMultistartInit(...).\nsolver = :rgd: Sets the algorithm for refinement. Possible options are:\n:rgd (default): Riemannian gradient descent.\n:als: Alternating least squares.\n:rcg: Riemannian conjugate gradient.\n:lbfgs: Limited-memory quasi-Newton refinement.\n\nExtended Options\n\ninit_point = nothing: Explicit initial point. If provided, it overrides the default initial point.\nwarm_init = BTDHOSVDMultistartInit(...): Searches for initial points using HOSVD for :alswarm, optionally screens them with short ALS runs, and returns the lowest-cost candidate.\nwarm_steps = 200: Once finding the best initial point, it runs this many ALS iterations to refine the initial point.\nwarm_block_method = :hooi or :sthosvd: Block update method used during warm start.\nwarm_block_maxiter = 20: Maximum number of inner iterations for each block update during warm start.\nwarm_rel_error_gate = 5e-2: Skips manifold refinement if the warm-start error is above this threshold.\nmaxiter = 500: Maximum number of Riemannian gradient descent iterations.\nstepsize = 0.01: Initial step size for line search in Riemannian gradient descent.\ntol = 1e-6: Convergence tolerance.\ngradient_mode = :riemannian: rgrad can be directly applied for manifold solvers. \nIf the model has a direct rgrad, it uses that.\nOtherwise it computes egrad and projects it to the tangent space.\nThis behavior is in src/solvers/abstract.jl (line 289).\nverbose = true: Enables progress output.\nblock_method = :hooi or :sthosvd: Block update method used for manifold solvers.\nblock_maxiter = 30: Maximum number of inner block-update iterations for manifold solvers.\nbtd_als_polish_maxiter = nothing: Number of final ALS polishing iterations. If nothing, an automatic budget is selected.\nThese settings are a robustness/quality feature for BTD-ALS. They are not used for manifold solvers.\nmax_stagnation_restarts = 1: More retries after a bad stagnated ALS pass. Higher values can improve solution quality, but increases runtime.\nstagnation_rel_error = 1e-4: This is a “bad final error” cutoff, not an improvement threshold. Lower value means restarts trigger more easily.\nrestart_candidates = 24: It usually improves chance of finding a better basin, but cost grows roughly linearly.\nrestart_screening_steps = 5: It uses this many quick ALS steps to screen restart candidates.\nrestart_block_maxiter = 20: Inner block-update limit during restart screening for BTD-ALS.\nrestart_seed = nothing: Optional random seed for restart generation for BTD-ALS.\n\nExample\n\njulia> using Random\njulia> Random.seed!(0)\njulia> A = randn(20, 15, 10); blocks = 10; ranks = (5, 4, 3)\njulia> res = btd(A, blocks, ranks; verbose = false)\nBTDResult{Float64}\n Blocks: 10\n Rel. error: 0.2625821087015455\n\n\n\n\n\n","title":"TensorKitchen.btd"},{"category":"type","location":"btd/#TensorKitchen.BTDResult","page":"BTD","text":"BTDResult{T}\n\nResult of block-term decomposition (btd); block components expose Tucker structure through accessors like core(blk), factors(blk), and blk.tensor.\n\nsolver_info: solver-specific diagnostics/metadata (NamedTuple). Typical keys include BTD-ALS restart diagnostics (total_iterations, stagnation_restarts, restart_rel_error_history) and BTD-TSD run settings (schedule, block_repeats, block_count, stepsize).\n\n\n\n\n\n","title":"TensorKitchen.BTDResult"},{"category":"function","location":"btd/#TensorKitchen.blocks","page":"BTD","text":"blocks(r::BTDResult)\n\nReturn the Tucker block components of a block-term decomposition result.\n\n\n\n\n\n","title":"TensorKitchen.blocks"},{"category":"method","location":"btd/#TensorKitchen.reconstruct-Tuple{BTDResult}","page":"BTD","text":"reconstruct(res::BTDResult)\n\nReconstruct the dense tensor represented by a block-term decomposition result by summing the reconstructed Tucker blocks.\n\n\n\n\n\n","title":"TensorKitchen.reconstruct"},{"category":"section","location":"tucker/#Tucker-Decomposition","page":"Tucker","text":"Approximating A by a Tucker decomposition\n\nhat A = C times_1 U times_2 V times_3 W\n\nwith multilinear rank mlrank can be computed as follows.\n\njulia> mlrank = (5, 4, 3)\njulia> tucker_res = tucker(A, mlrank)\nTuckerResult{Float64, 3}\n Original size: (20, 15, 10)\n Core size: (5, 4, 3)\n Multilinear rank: (5, 4, 3)\n Compression: 12.0x\n\nThe core C and the factor matrices (U V W) of the decomposition can be accessed as follows.\n\ncore(tucker_res)\nfactors(tucker_res)","title":"Tucker Decomposition"},{"category":"section","location":"tucker/#Tucker-Docs","page":"Tucker","text":"","title":"Tucker Docs"},{"category":"function","location":"tucker/#TensorKitchen.tucker","page":"Tucker","text":"tucker(A, ranks; method = :sthosvd, kwargs...) returns a TuckerResult\n\nComputes a Tucker decomposition of A with multilinear rank ranks.\n\nMain Options\n\nmethod = :sthosvd: Sets the Tucker decomposition algorithm. Possible options are:\n:sthosvd (default): Sequentially Truncated HOSVD. A direct one-pass decomposition, mainly used as a fast standalone Tucker approximation or as the default initializer for :hooiFast, deterministic, and usually a good initial point.\n:hooi: High-Order Orthogonal Iteration. Iteratively refines the Tucker factors, initialized by ST-HOSVD by default.\n\nExtended Options\n\nFor method = :hooi:\n\nmaxiter = 50: Maximum number of HOOI iterations.\ntol = 1e-8: Convergence tolerance based on change in relative reconstruction error.\ninit = :sthosvd\n:sthosvd: Uses ST-HOSVD to initialize the Tucker factors.\nTuckerResult: Uses an existing Tucker decomposition as the initial point.\n\nExample\n\njulia> using Random\njulia> Random.seed!(0)\njulia> A = randn(20, 15, 10); ranks = (5, 4, 3)\njulia> res = tucker(A, ranks; verbose = false)\nTuckerResult{Float64, 3}\n Original size: (20, 15, 10)\n Core size: (5, 4, 3)\n Multilinear rank: (5, 4, 3)\n Compression: 12.0x\n\n\nThe core tensor and factor matrices can be accessed by\n\ncore(res)\nfactors(res)\n\nA tensor approximation can be reconstructed by\n\nreconstruct(res)\n\nor equivalently\n\nreconstruct_tucker(core(res), factors(res))\n\nNotes\n\n:sthosvd is not an iterative solver. It directly returns a TuckerResult and does not expose solver-style outputs such as iteration counts or convergence diagnostics.\n:hooi is the iterative refinement method in the current Tucker implementation.\nImportant distinction: For the current Tucker implementation, do not use solver = :rgd, init = :auto, or manifold solvers like RGD, RCG, LBFGS, etc.\n\n\n\n\n\n","title":"TensorKitchen.tucker"},{"category":"type","location":"tucker/#TensorKitchen.TuckerResult","page":"Tucker","text":"TuckerResult{T, N}\n\nStores a Tucker decomposition: core tensor and factor matrices.\n\ncore::Array{T,N} — core tensor\nfactors::Vector{Matrix{T}} — orthonormal factor matrices\nprocessing_order::Vector{Int} — order modes were processed\nsingular_values::Vector{Vector{T}} — singular values per truncation\n\n\n\n\n\n","title":"TensorKitchen.TuckerResult"},{"category":"method","location":"tucker/#TensorKitchen.core-Tuple{TuckerResult}","page":"Tucker","text":"core(td::TuckerResult)\n\nReturn the Tucker core tensor for a Tucker result.\n\n\n\n\n\n","title":"TensorKitchen.core"},{"category":"method","location":"tucker/#TensorKitchen.factors-Tuple{TuckerResult}","page":"Tucker","text":"factors(td::TuckerResult)\n\nReturn the Tucker factor matrices.\n\n\n\n\n\n","title":"TensorKitchen.factors"},{"category":"method","location":"tucker/#TensorKitchen.multilinear_rank-Tuple{TuckerResult}","page":"Tucker","text":"multilinear_rank(td::TuckerResult)\n\nReturn the Tucker multilinear rank tuple, i.e. the size of the core tensor.\n\n\n\n\n\n","title":"TensorKitchen.multilinear_rank"},{"category":"method","location":"tucker/#TensorKitchen.factor_dims-Tuple{TuckerResult}","page":"Tucker","text":"factor_dims(td::TuckerResult)\n\nReturn the original mode dimensions represented by the Tucker factor matrices.\n\n\n\n\n\n","title":"TensorKitchen.factor_dims"},{"category":"method","location":"tucker/#TensorKitchen.reconstruct-Tuple{TuckerResult}","page":"Tucker","text":"reconstruct(td::TuckerResult)\n\nReconstruct the dense tensor represented by a Tucker decomposition result.\n\nFor a Tucker result with core S and factors U_1, ..., U_d, this returns S ×_1 U_1 ×_2 U_2 ... ×_d U_d.\n\n\n\n\n\n","title":"TensorKitchen.reconstruct"},{"category":"section","location":"cpd/#CPD","page":"CPD","text":"Here is how to approximate a tensor A by a CPD of rank r.\n\njulia> using TensorKitchen\njulia> A = randn(20, 15, 10)\njulia> r = 35\njulia> res = cpd(A, r)\nCPDResult{Float64}\n Order: 3\n Dimensions: (20, 15, 10)\n Rank: 35\n Rel. error: 0.4359141301703327\n\nNow, res contains a CP approximation of the 3-way tensor A,\n\nhat A = sum_i=1^r lambda_i a_i otimes b_i otimes c_i\n\nIt approximates A with relative error about 0.436.\n\nWe access the decomposition as follows.\n\nλ = weights(res)\nU = factors(res)\n\nHere, U is a triple of matrices (ABC), where the columns of A are the a_i and so on. These are called factor matrices.\n\nWe get the whole reconstructed tensor by \n\n = reconstruct(res)","title":"CPD"},{"category":"section","location":"cpd/#CPD-Docs","page":"CPD","text":"","title":"CPD Docs"},{"category":"function","location":"cpd/#TensorKitchen.cpd","page":"CPD","text":"cpd(A, r; kwargs...)\n\nComputes a rank-r CP approximation of A in two steps: (1) the first step finds an initial point; (2) the second step refines the initial point. Returns a CPDResult. If r is omitted, uses the smallest tensor mode as a heuristic rank.\n\nMain Options\n\ninit = :auto: Sets the algorithm to find the initial point. Possible options are:\n:auto: Uses a default CPD initializer. For solver = :als, this uses TuckerInit; otherwise, it uses an ALS warm start.\n:alswarm: Runs ALS first and uses the result as the initial point for refinement.\ncustomized initial point:\n:tucker (default when solver = :als): Uses a default Tucker initializer.\n:random: Uses a random initial point.\n:hosvd: Uses a HOSVD initial point.\nsolver = :rgd: Sets the algorithm for refinement. Possible options are:\nrgd (default): Riemannian gradient descent\nrgd_fixed: Riemannian gradient descent with fixed step size\nrcg: Riemannian conjugate gradient\nals: Alternating Least Squares\n\nExtended Options\n\np0 = nothing: Explicit initial point. If provided, it overrides the default initial point.\n:alswarm: ALS warm start option.\nwarm_init = TuckerInit(): Before finding the warm start initial point, this sets the good starting point for ALS.\nwarm_steps = 500: Once finding the best initial point from warm_init, it runs this many ALS iterations to refine the initial point.\nmaxiter = 500: Maximum number of Riemannian gradient descent iterations.\nstepsize = 1.0: Initial step size for line search in Riemannian gradient descent.\ntol = 1e-6: Convergence tolerance.\ngradient_mode = :riemannian: Gradient rule for manifold solvers. \nIf the model has a direct rgrad, it uses that.\nOtherwise it computes egrad and projects it to the tangent space.\nThis behavior is in src/solvers/abstract.jl (line 289).\ngeometry = :canonical: Sets the geometry of the manifold. Possible options are:\n:canonical: Standard CPD parameterization with the usual Euclidean factors and canonical Riemannian gradient handling. Best default for general unconstrained CPD.\n:squaring_metric: Nonnegative geometry based on squared latent coordinates. Enforces nonnegativity indirectly, but can become ill-conditioned near zero.\n:softplus_metric: Nonnegative geometry uses a regularized pullback-inspired geometry induced by the softplus chart. Smoother and usually more stable near zero than :squaring_metric.\n:native: Native CP manifold geometry using the model’s intrinsic CP/Segre representation not for nonnegative=true. Best for structured join layouts with Manifolds.Segre summands.\nverbose = true: Enables progress output.\nnonnegative::Bool = false: Nonnegative CPD option to be selected by the user. (same as nncpd)\npullback_eps = 1e-8: Regularization parameter for pullback-style nonnegative geometries.\n\nNotes\n\nsolver = :als does not use manifold geometry. In that case:\ngeometry must be :canonical\ngradient_mode is ignored except for validation\n:squaring_metric and :softplus_metric require nonnegative = true.\nWhen nonnegative = true, cpd(...) routes to nncpd(...). In that route:\nif solver != :als and geometry is left at :canonical, the effective geometry becomes :softplus_metric\nif stepsize is left at 1.0, the effective default becomes 0.01\nif init = :tucker, the effective initializer becomes :alswarm\n\nExample\n\njulia> A = randn(20, 15, 10); r = 35\njulia> res = cpd(A, r)\nCPDResult{Float64}\n Order: 3\n Dimensions: (20, 15, 10)\n Rank: 35\n Rel. error: 0.4359141301703327\n\n\n\n\n\n","title":"TensorKitchen.cpd"},{"category":"function","location":"cpd/#TensorKitchen.nncpd","page":"CPD","text":" nncpd(A, r; kwargs...)\n\nComputes a nonnegative rank-r CP approximation of A in two steps: (1) the first step finds an initial point; (2) the second step refines the initial point. Returns a CPDResult. If r is omitted, uses the smallest tensor mode as a heuristic rank. cpd(A, r; nonnegative=true, ...) routes here and adopts the same effective defaults.\n\nOptions\n\nThe options are the same as for cpd.\n\nGeometry guide:\n\ngeometry=:softplus_metric Default and usually the safest choice. \ngeometry=:squaring_metric Uses a regularized pullback-inspired geometry induced by the squaring chart.\ngeometry=:canonical Plain nonnegative CP coordinates without the pullback-style manifold geometry. This is the natural choice with solver=:als.\n\nExample\n\njulia> A = randn(20, 15, 10); r = 35;\njulia> B = abs.(A)\njulia> nncpd(B, r)\nCPDResult{Float64}\n Order: 3\n Dimensions: (20, 15, 10)\n Rank: 35\n Rel. error: 0.3765605093526155\n\n\n\n\n\n","title":"TensorKitchen.nncpd"},{"category":"type","location":"cpd/#TensorKitchen.CPDResult","page":"CPD","text":"CPDResult{T}\n\nResult of a Canonical Polyadic Decomposition.\n\nStores the decoded CP representation together with solver diagnostics:\n\ncomponents: rank-one tensor components\nweights: component weights\nfactors: factor matrices\ncost: final objective function value at the returned solution\nrel_error: final relative reconstruction error\ngrad_norm: norm of the final optimization gradient reported by the solver; for manifold solvers this is the Riemannian gradient norm\niterations: number of refinement iterations\nconverged: whether the solver reported convergence\nsolver: solver optimization method used to produce the result\nsolver_info: solver-specific diagnostics/metadata (NamedTuple). Typical keys include:\ninitial_stepsize_eff (RGD), memory_size (LBFGS), cautious_update (LBFGS),\ninitial_scale, linesearch, has_preconditioner (LBFGS), and nncp_pullback_eps (NNCP).\n\n\n\n\n\n","title":"TensorKitchen.CPDResult"},{"category":"method","location":"cpd/#TensorKitchen.weights-Tuple{CPDResult}","page":"CPD","text":"weights(r::CPDResult)\n\nReturn the CP component weights stored in a CPD result.\n\n\n\n\n\n","title":"TensorKitchen.weights"},{"category":"method","location":"cpd/#TensorKitchen.factors-Tuple{CPDResult}","page":"CPD","text":"factors(res::CPDResult)\n\nReturn the CP factor matrices of res as a vector [U₁, U₂, ..., U_N], where each U_m has size size(A, m) × rank.\n\n\n\n\n\n","title":"TensorKitchen.factors"},{"category":"method","location":"cpd/#TensorKitchen.reconstruct-Tuple{CPDResult}","page":"CPD","text":"reconstruct(res::CPDResult)\n\nReconstruct the dense tensor represented by a CP decomposition result.\n\nFor a rank-R CPD result, this returns sum(weights(res)[k] * u_1k ⊗ ... ⊗ u_Nk for k = 1:R).\n\n\n\n\n\n","title":"TensorKitchen.reconstruct"},{"category":"section","location":"#TensorKitchen.jl:-tensor-decompositions-in-Julia","page":"Home","text":"\n\nTensorKitchen.jl is a Julia package for tensor decompositions.\n\nCPD documentation\nTucker documentation\nBTD documentation\nJoin decomposition documentation\nUtilities\nPipeline\nReferences","title":"TensorKitchen.jl: tensor decompositions in Julia"},{"category":"section","location":"#Notes","page":"Home","text":"The package is currently an early version and will be updated frequently in the near future.\n\nThe implementation is based on combining algebraic algorithms like ALS (see, e.g., the textbook by Kolda and Ballard) and Riemannian optimization from Manopt.jl.\n\nWhat currently works is \n\nCanonical Polyadic Decomposition (CPD)\nTucker Decomposition\nNonnegative Canonical Polyadic Decomposition (NNCPD)\nBlock Term Decomposition (BTD)\nJoin Decompositions\n\n\n\nThe next updates will include \n\nHandling of swamps/plateaus in the optimization step\nDocumentation\nImproved User Interface\nGPU Support \nLL1 Decomposition (3-way specialized BTD)\nSymmetric CP / Waring Decomposition\nPartially Symmetric CP\nTensor Trains","title":"Notes"},{"category":"section","location":"utils/#General-Utilities","page":"Utilities","text":"","title":"General Utilities"},{"category":"function","location":"utils/#TensorKitchen.save_result","page":"Utilities","text":"save_result(path::AbstractString, result)\n\nSave a TensorKitchen result object to path for later use.\n\nThis uses Julia's built-in Serialization format, so the file is a Julia-native binary file. \nIt is intended for saving results during local experiments, benchmarks, and development workflows.\nCommon file extensions are .jls, .julia, .bin, or .tkresult; the extension is only a convention.\n\nExamples\n\nA = randn(20, 15, 10)\nres = cpd(A, 35)\n\nsave_result(\"cpd_rank35.jls\", res)\n\nYou can also save richer records by serializing a NamedTuple of the result and additional metadata.\n\nrecord = (\n method = :cpd,\n rank = 35,\n input_size = size(A),\n result = res,\n)\nsave_result(\"experiment_cpd_rank35.jls\", record)\n\n\n\n\n\n","title":"TensorKitchen.save_result"},{"category":"function","location":"utils/#TensorKitchen.load_result","page":"Utilities","text":"load_result(path::AbstractString)\n\nLoad a previously saved TensorKitchen result or experiment record from path for later use.\n\nThe file must be written with save_result, or otherwise created using Julia's Serialization.serialize.\n\nExamples\n\nres = load_result(\"cpd_rank35.jls\")\n\nweights(res)\nfactors(res)\nreconstruct(res)\n\nIf the saved object was an experiment record, you can access the stored fields:\n\nrecord = load_result(\"experiment_cpd_rank35.jls\")\n\nrecord.method\nrecord.rank\nrecord.input_size\nrecord.result\nrecord.result.weights\nrecord.result.factors\nreconstruct(record.result)\n\n\n\n\n\n","title":"TensorKitchen.load_result"}] } diff --git a/docs/src/PIPELINE.md b/docs/src/PIPELINE.md index 7b1f74a..b1ef8bd 100644 --- a/docs/src/PIPELINE.md +++ b/docs/src/PIPELINE.md @@ -5,11 +5,11 @@ converters. ## Public entry points -- `cpd(A, r; ...)` -> `CPDResult` -- `nncpd(A, r; ...)` -> `CPDResult` -- `btd(A, blocks, ranks; ...)` -> `BTDResult` -- `tucker(A, ranks; method=...)` -> `TuckerResult` -- `approx(...)` -> `ApproxResult` or auto-routed `CPDResult`/`BTDResult` +- CP Decomposition `cpd(A, r; ...)` +- Nonnegative CP Decomposition `nncpd(A, r; ...)` +- Block Term Decomposition `btd(A, blocks, ranks; ...)` +- Tucker Decomposition `tucker(A, ranks; method=...)` +- Join Decomposition `approx(...)` ## Default behavior (quick reference) @@ -30,8 +30,7 @@ converters. - `tucker(A, ranks)`: - `method = :sthosvd` - `approx(model::JoinModel)`: - - `init = :alswarm` - - `warm_steps = 500` + - `init = :random` - `solver = :rgd` ## Core execution architecture @@ -54,13 +53,14 @@ symbol-to-solver dispatch layer (`:rgd`, `:rcg`, `:lbfgs`, `:als`, `:btd_tsd`). 1. Build `JoinModel(A, r; geometry=...)` with `CPDBackend` 2. Normalize/validate options (`solver`, `geometry`, `gradient_mode`, normalization policy) 3. Solve through `_solve_model(...)` -4. Optionally run nonnegative ALS polishing (for selected nonnegative paths) -5. Convert to `CPDResult` +4. Convert to `CPDResult` Notes: - `:als` means CP-ALS. - Manifold solvers (`:rgd`, `:rgd_fixed`, `:rcg`, `:lbfgs`) share dispatch with other pipelines. +- For `solver != :als`, `init = :auto` resolves to `:alswarm`, so CPD and NNCPD start from an ALS warm point before manifold refinement. +- Generic `approx(...)` does not use CPD's ALS warm-start path unless it auto-routes to `cpd(...)`. ### BTD (`btd`) @@ -76,7 +76,7 @@ Notes: 6. If `solver != :als`, optionally polish with BTD-ALS (`btd_als_polish_maxiter`) 7. Convert to `BTDResult` -Polish step usefulness (brief): +Polish step usefulness: - Usually helpful for a small final `rel_error` reduction after RGD converges near a good basin. - Most useful for quality-focused runs (benchmarks, final fits). @@ -93,7 +93,7 @@ BTD-ALS stabilization behavior: - Tracks per-iteration fit change (`|rel_t - rel_{t-1}|`) - Detects stagnation when fit change is tiny but `rel_error` remains high - Can restart from fresh multistart pool (`max_stagnation_restarts`) -- Reports true final Riemannian gradient norm (`grad_norm`) instead of a placeholder +- Reports true final Riemannian gradient norm (`grad_norm`) ### Tucker (`tucker`) @@ -113,6 +113,11 @@ It dispatches directly to decomposition routines: `dispatch=:cpd`, `:btd`, and `:generic` force behavior. +For the generic `JoinModel` path, `approx(...)` starts from `init = :random` +by default and then runs the selected manifold solver. It does not run an ALS +warm-start stage, because a general join component does not necessarily expose +factor matrices or least-squares block updates. + ## Result types and post-processing - `CPDResult` @@ -127,7 +132,4 @@ Common utilities: ## File map -- API entry points: `src/api/approx.jl`, `src/api/cpd.jl`, `src/api/nncpd.jl`, `src/api/btd.jl` -- Routing helpers: `src/dispatch/approx_routing.jl` -- Solver dispatch core: `src/solvers/solve_dispatch.jl` -- BTD backend/init details: `src/btd/model.jl`, `src/solvers/btd_als.jl` +- API entry points: `src/api/approx.jl`, `src/api/cpd.jl`, `src/api/nncpd.jl`, `src/api/btd.jl`, `src/api/tucker.jl` diff --git a/docs/src/index.md b/docs/src/index.md index 5404928..e0749f2 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -1,5 +1,6 @@ -# TensorKitchen.jl Documentation +# TensorKitchen.jl: tensor decompositions in Julia + **TensorKitchen.jl** is a Julia package for tensor decompositions. @@ -13,7 +14,7 @@ ## Notes -The package is currently at a pre-alpha stage. +The package is currently an early version and will be updated frequently in the near future. The implementation is based on combining algebraic algorithms like ALS (see, e.g., the [textbook by Kolda and Ballard](https://users.wfu.edu/ballard/pdfs/tensor_textbook.pdf)) and Riemannian optimization from [Manopt.jl](https://manoptjl.org/stable/). diff --git a/docs/tucker/index.html b/docs/tucker/index.html index cd7b853..8c905a0 100644 --- a/docs/tucker/index.html +++ b/docs/tucker/index.html @@ -1,5 +1,5 @@ -Tucker · Documentation
GitHub

Tucker Decomposition

Approximating A by a Tucker decomposition

\[\hat A = C \times_1 U \times_2 V \times_3 W\]

with multilinear rank mlrank can be computed as follows.

julia> mlrank = (5, 4, 3)
+Tucker · TensorKitchen.jl
GitHub
Tucker Decomposition
Approximating A by a Tucker decomposition
\[\hat A = C \times_1 U \times_2 V \times_3 W\]
with multilinear rank mlrank can be computed as follows.
julia> mlrank = (5, 4, 3)
 julia> tucker_res = tucker(A, mlrank)
 TuckerResult{Float64, 3}
   Original size:    (20, 15, 10)
@@ -16,6 +16,4 @@
   Multilinear rank: (5, 4, 3)
   Compression:      12.0x
 
The core tensor and factor matrices can be accessed by
core(res)
-factors(res)
A tensor approximation can be reconstructed by
reconstruct(res)
or equivalently
reconstruct_tucker(core(res), factors(res))
Notes
  • :sthosvd is not an iterative solver. It directly returns a TuckerResult and does not expose solver-style outputs such as iteration counts or convergence diagnostics.
  • :hooi is the iterative refinement method in the current Tucker implementation.
  • Important distinction: For the current Tucker implementation, do not use solver = :rgd, init = :auto, or manifold solvers like RGD, RCG, LBFGS, etc.
source
TensorKitchen.TuckerResult — Type
TuckerResult{T, N}

Stores a Tucker decomposition: core tensor and factor matrices.

  • core::Array{T,N} — core tensor
  • factors::Vector{Matrix{T}} — orthonormal factor matrices
  • processing_order::Vector{Int} — order modes were processed
  • singular_values::Vector{Vector{T}} — singular values per truncation
source
TensorKitchen.reconstruct — Method
reconstruct(td::TuckerResult) reconstructs the tensor from Tucker decomposition
-
-A = S ×₁ U₁ ×₂ U₂ ⋯ ×_d U_d
source
+factors(res)

A tensor approximation can be reconstructed by

reconstruct(res)

or equivalently

reconstruct_tucker(core(res), factors(res))

Notes

source
TensorKitchen.TuckerResult — Type
TuckerResult{T, N}

Stores a Tucker decomposition: core tensor and factor matrices.

  • core::Array{T,N} — core tensor
  • factors::Vector{Matrix{T}} — orthonormal factor matrices
  • processing_order::Vector{Int} — order modes were processed
  • singular_values::Vector{Vector{T}} — singular values per truncation
source
TensorKitchen.core — Method
core(td::TuckerResult)

Return the Tucker core tensor for a Tucker result.

source
TensorKitchen.factors — Method
factors(td::TuckerResult)

Return the Tucker factor matrices.

source
TensorKitchen.multilinear_rank — Method
multilinear_rank(td::TuckerResult)

Return the Tucker multilinear rank tuple, i.e. the size of the core tensor.

source
TensorKitchen.factor_dims — Method
factor_dims(td::TuckerResult)

Return the original mode dimensions represented by the Tucker factor matrices.

source
TensorKitchen.reconstruct — Method
reconstruct(td::TuckerResult)

Reconstruct the dense tensor represented by a Tucker decomposition result.

For a Tucker result with core S and factors U_1, ..., U_d, this returns S ×_1 U_1 ×_2 U_2 ... ×_d U_d.

source
diff --git a/docs/utils/index.html b/docs/utils/index.html index 9562261..47acc80 100644 --- a/docs/utils/index.html +++ b/docs/utils/index.html @@ -1,2 +1,23 @@ -Utilities · Documentation
GitHub
+Utilities · TensorKitchen.jl
GitHub

General Utilities

TensorKitchen.save_result — Function
save_result(path::AbstractString, result)

Save a TensorKitchen result object to path for later use.

  • This uses Julia's built-in Serialization format, so the file is a Julia-native binary file.
  • It is intended for saving results during local experiments, benchmarks, and development workflows.
  • Common file extensions are .jls, .julia, .bin, or .tkresult; the extension is only a convention.

Examples

A = randn(20, 15, 10)
+res = cpd(A, 35)
+
+save_result("cpd_rank35.jls", res)
  • You can also save richer records by serializing a NamedTuple of the result and additional metadata.
record = (
+    method = :cpd,
+    rank = 35,
+    input_size = size(A),
+    result = res,
+)
+save_result("experiment_cpd_rank35.jls", record)
source
TensorKitchen.load_result — Function
load_result(path::AbstractString)

Load a previously saved TensorKitchen result or experiment record from path for later use.

  • The file must be written with save_result, or otherwise created using Julia's Serialization.serialize.

Examples

res = load_result("cpd_rank35.jls")
+
+weights(res)
+factors(res)
+reconstruct(res)

If the saved object was an experiment record, you can access the stored fields:

record = load_result("experiment_cpd_rank35.jls")
+
+record.method
+record.rank
+record.input_size
+record.result
+record.result.weights
+record.result.factors
+reconstruct(record.result)
source
diff --git a/src/api/cpd.jl b/src/api/cpd.jl index 119cd09..69c7701 100644 --- a/src/api/cpd.jl +++ b/src/api/cpd.jl @@ -399,9 +399,9 @@ If `r` is omitted, uses the smallest tensor mode as a heuristic rank. * `stepsize = 1.0`: Initial step size for line search in Riemannian gradient descent. * `tol = 1e-6`: Convergence tolerance. * `gradient_mode = :riemannian`: Gradient rule for manifold solvers. - - If the model has a direct rgrad, it uses that. - - Otherwise it computes egrad and projects it to the tangent space. - - This behavior is in src/solvers/abstract.jl (line 289). + - If the model has a direct rgrad, it uses that. + - Otherwise it computes egrad and projects it to the tangent space. + - This behavior is in src/solvers/abstract.jl (line 289). * `geometry = :canonical`: Sets the geometry of the manifold. Possible options are: - `:canonical`: Standard CPD parameterization with the usual Euclidean factors and canonical Riemannian gradient handling. Best default for general unconstrained CPD. - `:squaring_metric`: Nonnegative geometry based on squared latent coordinates. Enforces nonnegativity indirectly, but can become ill-conditioned near zero. diff --git a/src/core/tensor_ops.jl b/src/core/tensor_ops.jl index f999392..5f89be5 100644 --- a/src/core/tensor_ops.jl +++ b/src/core/tensor_ops.jl @@ -325,7 +325,7 @@ function rank1_mode_contract_column!( end return out end - +# 3D version function rank1_mode_contract_column( A::AbstractArray{T,3}, U::AbstractVector{<:AbstractMatrix{T}}, @@ -347,7 +347,7 @@ function rank1_mode_contract_column( end return out end - +# 4D version function rank1_mode_contract_column( A::AbstractArray{T,4}, U::AbstractVector{<:AbstractMatrix{T}}, @@ -483,7 +483,7 @@ function cp_inner_AX( ) end -@inline function _cp_residual_sq_from_gram_unreliable( +@inline function _cp_residual_sq_from_G_unreliable( n2::T, normA2::T, normX2::T, @@ -503,7 +503,7 @@ function cp_residual_stats( normX2 = cp_reconstruction_norm2(components) innerAX = cp_inner_AX(A, components) n2 = normA2 + normX2 - 2 * innerAX - if _cp_residual_sq_from_gram_unreliable(n2, normA2, normX2, innerAX) + if _cp_residual_sq_from_G_unreliable(n2, normA2, normX2, innerAX) return cp_residual_stats_explicit(A, normA2, components) end return (n2, T(0.5) * n2, _relative_error_frob_sq(n2, normA2)) diff --git a/src/cpd/core/mttkrp.jl b/src/cpd/core/mttkrp.jl index de993f0..47be706 100644 --- a/src/cpd/core/mttkrp.jl +++ b/src/cpd/core/mttkrp.jl @@ -1,8 +1,8 @@ # cpd/core/mttkrp.jl — CPD-specific MTTKRP kernels and dispatch -export mttkrp, khatri_rao -# --------------------------------------------------------------------------- -# Auto policy -# --------------------------------------------------------------------------- +# Improvement of resolving mttkrp bottleneck still in progress +export mttkrp, mttkrp!, khatri_rao, khatri_rao! +@inline _mttkrp_needs_kr_workspace(method::Symbol) = method == :khatri_rao +@inline _mttkrp_needs_tmp_workspace(method::Symbol) = method in (:direct3, :direct4) @inline function _mttkrp_auto_method_3way(kr_rows::Int, r::Int, mode::Int) # Benchmark-guided 3-way table: @@ -32,10 +32,28 @@ end end end -# --------------------------------------------------------------------------- -# Khatri-Rao helpers -# --------------------------------------------------------------------------- +@inline function _mttkrp_resolve_method( + method::Symbol, + dims::NTuple{N,Int}, + r::Int, + mode::Int, +) where {N} + method == :auto && return _mttkrp_auto_method(dims, r, mode) + method == :khatri_rao && return :khatri_rao + if method == :direct + N == 3 && return :direct3 + N == 4 && return :direct4 + return :contract + else + throw( + ArgumentError( + "Unknown mttkrp method=$method. Use :auto, :khatri_rao, or :direct.", + ), + ) + end +end +# forming Khatri-Rao product helper, the loop is costly. function khatri_rao(mats::AbstractVector{<:AbstractMatrix{T}}) where {T<:AbstractFloat} if isempty(mats) throw(ArgumentError("khatri_rao: empty matrix list")) @@ -142,10 +160,7 @@ function _mttkrp_khatri_rao!( return out end -# --------------------------------------------------------------------------- # Direct and contraction kernels -# --------------------------------------------------------------------------- - @inline function _accumulate_scaled_columns!( out::AbstractMatrix{T}, tmp::AbstractMatrix{T}, @@ -324,10 +339,7 @@ function _mttkrp_contract( return out end -# --------------------------------------------------------------------------- -# Public API -# --------------------------------------------------------------------------- - +#Public API function mttkrp( A::AbstractArray{T,N}, components::Vector{RankOneTensor{T}}, @@ -361,7 +373,7 @@ function mttkrp( throw(DimensionMismatch("mttkrp: all factors must have same column count")) end - method_eff = method == :auto ? _mttkrp_auto_method(dims, r, mode) : method + method_eff = _mttkrp_resolve_method(method, dims, r, mode) if method_eff == :khatri_rao return _mttkrp_khatri_rao(A, U, mode) @@ -378,7 +390,7 @@ function mttkrp( else throw( ArgumentError( - "Unknown mttkrp method=$method. Use :auto, :khatri_rao, :direct3, :direct4, or :contract.", + "Unknown mttkrp method=$method. Use :auto, :khatri_rao, or :direct.", ), ) end @@ -395,14 +407,30 @@ function mttkrp!( kr_work = nothing, ) where {T<:AbstractFloat,N} dims = size(A) + mode < 1 && throw(ArgumentError("mode must be >= 1")) + mode > N && throw(ArgumentError("mode must be <= ndims(A)")) + isempty(U) && throw(ArgumentError("mttkrp: factor list is empty")) + length(U) == N || + throw(DimensionMismatch("mttkrp: expected $N factor matrices, got $(length(U))")) + + r = size(U[1], 2) + @inbounds for m = 1:N + size(U[m], 1) == dims[m] || throw( + DimensionMismatch( + "mttkrp: U[$m] has $(size(U[m], 1)) rows, expected $(dims[m])", + ), + ) + size(U[m], 2) == r || + throw(DimensionMismatch("mttkrp: all factors must have same column count")) + end + size(out, 1) == dims[mode] || throw( DimensionMismatch("mttkrp!: out has $(size(out,1)) rows, expected $(dims[mode])"), ) - r = size(U[1], 2) size(out, 2) == r || throw(DimensionMismatch("mttkrp!: out has $(size(out,2)) columns, expected $r")) - method_eff = method == :auto ? _mttkrp_auto_method(dims, r, mode) : method + method_eff = _mttkrp_resolve_method(method, dims, r, mode) if method_eff == :khatri_rao isnothing(kr_buf) && throw( ArgumentError("mttkrp!: method=:khatri_rao requires a KR workspace buffer"), @@ -440,7 +468,7 @@ function mttkrp!( else throw( ArgumentError( - "Unknown mttkrp method=$method. Use :auto, :khatri_rao, :direct3, :direct4, or :contract.", + "Unknown mttkrp method=$method. Use :auto, :khatri_rao, or :direct.", ), ) end diff --git a/src/results/reconstruct.jl b/src/results/reconstruct.jl index 3a60cf2..98da321 100644 --- a/src/results/reconstruct.jl +++ b/src/results/reconstruct.jl @@ -1,11 +1,21 @@ # results/reconstruct.jl — result reconstruction helpers export reconstruct +""" + reconstruct(res::CPDResult) + +Reconstruct the dense tensor represented by a CP decomposition result. + +For a rank-`R` CPD result, this returns +`sum(weights(res)[k] * u_1k ⊗ ... ⊗ u_Nk for k = 1:R)`. +""" reconstruct(res::CPDResult) = reconstruct_cpd_rankr(components(res)) """ reconstruct(res::ApproxResult) - reconstruct(res::BTDResult) + +Reconstruct the dense ambient object represented by a generic join +approximation result by summing its component tensors. """ function reconstruct(res::ApproxResult) comps = components(res) @@ -18,6 +28,12 @@ function reconstruct(res::ApproxResult) return X end +""" + reconstruct(res::BTDResult) + +Reconstruct the dense tensor represented by a block-term decomposition result +by summing the reconstructed Tucker blocks. +""" function reconstruct(res::BTDResult) comps = components(res) isempty(comps) && throw(ArgumentError("BTDResult has no components to reconstruct.")) diff --git a/src/solvers/btd_tsd.jl b/src/solvers/btd_tsd.jl index 4fc38b1..3f12788 100644 --- a/src/solvers/btd_tsd.jl +++ b/src/solvers/btd_tsd.jl @@ -1,4 +1,16 @@ # solvers/btd_tsd.jl — blockwise tangent-subspace descent for BTD + +# This solver is inspired by tangent-subspace descent (TSD) in +# Gutman and Ho-Nguyen, "Coordinate Descent Without Coordinates: +# Tangent Subspace Descent on Riemannian Manifolds", +# Mathematics of Operations Research, 48(1):127–159, 2023. +# +# In the present implementation, the selected tangent subspace is the +# tangent space of one Tucker block inside the BTD join/product model. +# Thus this is a BTD-specialized blockwise TSD / Riemannian block-coordinate +# descent method, rather than the exact algorithm studied for Stiefel/orthogonal +# examples in the original TSD paper. + export BTDTSDSolver, TSDSolver """ @@ -8,6 +20,11 @@ export BTDTSDSolver, TSDSolver Blockwise tangent-subspace descent specialized to BTD. Each block update uses the projected Tucker-block tangent direction and accepts it only when a block Armijo decrease condition is satisfied. + +This solver is inspired by tangent-subspace descent (TSD) in +Gutman and Ho-Nguyen, "Coordinate Descent Without Coordinates: +Tangent Subspace Descent on Riemannian Manifolds", +Mathematics of Operations Research, 48(1):127–159, 2023. """ struct BTDTSDSolver <: AbstractFirstOrderSolver stepsize::Float64 @@ -18,12 +35,27 @@ struct BTDTSDSolver <: AbstractFirstOrderSolver armijo_alpha_min::Float64 end -const TSDSolver = BTDTSDSolver +""" + BTDTSDSolver(; schedule=:cyclic, block_repeats=1, ...) + +Blockwise tangent-subspace descent for BTD. +The `schedule` keyword controls the order in which Tucker blocks are updated. + +- `schedule = :cyclic`: updates blocks in the deterministic order + `1, 2, ..., R` at every sweep. + +- `schedule = :random`: updates all blocks once per sweep, but in a fresh + random permutation. This is random reshuffling, not sampling with replacement. + +The cyclic schedule is the default because it is deterministic and easiest to +analyze. The random schedule can be useful experimentally when different BTD +blocks compete strongly for the same residual structure. +""" function BTDTSDSolver(; stepsize::Real = 1.0, schedule::Symbol = :cyclic, - block_repeats::Int = 1, + block_repeats::Int = 1, # number of times to repeat the block update sequence at each iteration armijo_contraction::Real = 0.5, armijo_sufficient_decrease::Real = 1e-4, armijo_alpha_min::Real = 1e-12, diff --git a/src/solvers/cp_als.jl b/src/solvers/cp_als.jl index b908dea..3135f84 100644 --- a/src/solvers/cp_als.jl +++ b/src/solvers/cp_als.jl @@ -3,13 +3,13 @@ export ALSSolver, fit_cp_als struct CPALSWorkspace - grams::Any - V::Any + Gs::Any # Unpacked G matrices (U[n] * U[n]') + V::Any # V matrix (U[1] * U[1]') - cached MTTKRP result transposed_work::Any - denom_work::Any + denom_work::Any # Denominator matrix for non-LS updates mttkrp_bufs::Any mttkrp_tmp_work::Any - mttkrp_kr_work::Any + mttkrp_kr_work::Any # MTTKRP kernel result buffer mttkrp_kr_work2::Any cross_buf::Any end @@ -59,20 +59,34 @@ end function CPALSWorkspace( A::AbstractArray{T}, dims::NTuple{N,Int}, - r::Int, + r::Int; + mttkrp_method::Symbol = :auto, ) where {T<:AbstractFloat,N} - grams = [_cp_als_matrix_workspace_like(A, r, r) for _ = 1:N] + Gs = [_cp_als_matrix_workspace_like(A, r, r) for _ = 1:N] V = _cp_als_matrix_workspace_like(A, r, r) transposed_work = [_cp_als_matrix_workspace_like(A, r, dims[n]) for n = 1:N] denom_work = [_cp_als_matrix_workspace_like(A, dims[n], r) for n = 1:N] mttkrp_bufs = [_cp_als_matrix_workspace_like(A, dims[n], r) for n = 1:N] - mttkrp_tmp_work = [_cp_als_matrix_workspace_like(A, dims[n], r) for n = 1:N] - kr_rows = [div(prod(dims), dims[n]) for n = 1:N] - mttkrp_kr_work = [_cp_als_matrix_workspace_like(A, kr_rows[n], r) for n = 1:N] - mttkrp_kr_work2 = [_cp_als_matrix_workspace_like(A, kr_rows[n], r) for n = 1:N] + total_dim_prod = prod(dims) + resolved_mttkrp_methods = + [_mttkrp_resolve_method(mttkrp_method, dims, r, n) for n = 1:N] + mttkrp_tmp_work = Any[ + _mttkrp_needs_tmp_workspace(resolved_mttkrp_methods[n]) ? + _cp_als_matrix_workspace_like(A, dims[n], r) : nothing for n = 1:N + ] + mttkrp_kr_work = Any[ + _mttkrp_needs_kr_workspace(resolved_mttkrp_methods[n]) ? + _cp_als_matrix_workspace_like(A, div(total_dim_prod, dims[n]), r) : nothing for + n = 1:N + ] + mttkrp_kr_work2 = Any[ + _mttkrp_needs_kr_workspace(resolved_mttkrp_methods[n]) ? + _cp_als_matrix_workspace_like(A, div(total_dim_prod, dims[n]), r) : nothing for + n = 1:N + ] cross_buf = _cp_als_matrix_workspace_like(A, r, r) return CPALSWorkspace( - grams, + Gs, V, transposed_work, denom_work, @@ -85,7 +99,7 @@ function CPALSWorkspace( end -@inline function _update_gram!( +@inline function _update_G!( G::AbstractMatrix{T}, U::AbstractMatrix{T}, ) where {T<:AbstractFloat} @@ -93,15 +107,15 @@ end return G end -@inline function _hadamard_gram_except!( +@inline function _hadamard_G_except!( V::AbstractMatrix{T}, - grams::AbstractVector{<:AbstractMatrix{T}}, + Gs::AbstractVector{<:AbstractMatrix{T}}, skip::Int, ) where {T<:AbstractFloat} fill!(V, one(T)) - @inbounds for m in eachindex(grams) + @inbounds for m in eachindex(Gs) m == skip && continue - V .*= grams[m] + V .*= Gs[m] end @inbounds for i in axes(V, 1) V[i, i] += eps(T) @@ -111,11 +125,11 @@ end function _solve_right_spd!( Udest::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, V::AbstractMatrix{T}, workT::AbstractMatrix{T}, ) where {T<:AbstractFloat} - copyto!(workT, transpose(G)) + copyto!(workT, transpose(M_mttkrp)) F = cholesky!(Hermitian(V)) ldiv!(F, workT) copyto!(Udest, transpose(workT)) @@ -124,24 +138,24 @@ end @inline function _als_mode_update!( Udest::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, V::AbstractMatrix{T}, workT::AbstractMatrix{T}, ) where {T<:AbstractFloat} - return _solve_right_spd!(Udest, G, V, workT) + return _solve_right_spd!(Udest, M_mttkrp, V, workT) end """ - _cp_als_normalize_weights_and_grams!(λ, U, grams, normalization_policy, update_policy) + _cp_als_normalize_weights_and_Gs!(λ, U, Gs, normalization_policy, update_policy) Post-sweep step: reset `λ` to one, apply the normalization policy, clamp to -nonnegative if required, and refresh Gram matrices. Single pass — see the +nonnegative if required, and refresh G matrices. Single pass — see the module-level design note. """ -@inline function _cp_als_normalize_weights_and_grams!( +@inline function _cp_als_normalize_weights_and_Gs!( λ, U, - grams, + Gs, normalization_policy, update_policy, ) @@ -149,8 +163,8 @@ module-level design note. fill!(λ, one(T)) normalize_components!(U, λ, normalization_policy) update_policy != :ls && _clamp_nonnegative!(λ, U) - @inbounds for n in eachindex(grams) - _update_gram!(grams[n], U[n]) + @inbounds for n in eachindex(Gs) + _update_G!(Gs[n], U[n]) end return nothing end @@ -160,7 +174,7 @@ function _cp_als_stats( normA2::T, λ::AbstractVector{T}, U::AbstractVector{<:AbstractMatrix{T}}, - grams::AbstractVector{<:AbstractMatrix{T}}; + Gs::AbstractVector{<:AbstractMatrix{T}}; mttkrp_method::Symbol = :auto, mttkrp_buf::Union{Nothing,AbstractMatrix{T}} = nothing, mttkrp_work::Union{Nothing,AbstractMatrix{T}} = nothing, @@ -169,14 +183,14 @@ function _cp_als_stats( cross_buf::Union{Nothing,AbstractMatrix{T}} = nothing, ) where {T<:AbstractFloat,N} if isnothing(cross_buf) - cross = copy(grams[1]) + cross = copy(Gs[1]) else cross = cross_buf - copyto!(cross, grams[1]) + copyto!(cross, Gs[1]) end - @inbounds for m = 2:length(grams) - cross .*= grams[m] + @inbounds for m = 2:length(Gs) + cross .*= Gs[m] end normX2 = zero(T) @inbounds for j in axes(cross, 2) @@ -205,7 +219,7 @@ function _cp_als_stats( innerAX += λ[k] * dot(@view(U[1][:, k]), @view(M1[:, k])) end n2 = normA2 + normX2 - 2 * innerAX - if _cp_residual_sq_from_gram_unreliable(n2, normA2, normX2, innerAX) + if _cp_residual_sq_from_G_unreliable(n2, normA2, normX2, innerAX) return cp_residual_stats_explicit(A, normA2, λ, U) end return (n2, T(0.5) * n2, _relative_error_frob_sq(n2, normA2)) @@ -271,10 +285,10 @@ function fit_cp_als( nnls_row_tol = sqrt(eps(T)) nnls_row_tol_min = nnls_row_tol / 10 - workspace = CPALSWorkspace(A, dims, r) - grams = workspace.grams + workspace = CPALSWorkspace(A, dims, r; mttkrp_method = mttkrp_method) + Gs = workspace.Gs @inbounds for n = 1:N - _update_gram!(grams[n], U[n]) + _update_G!(Gs[n], U[n]) end V = workspace.V transposed_work = workspace.transposed_work @@ -294,8 +308,8 @@ function fit_cp_als( pg_sq = zero(T) u_sq = zero(T) for n = 1:N - _hadamard_gram_except!(V, grams, n) - G = mttkrp!( + _hadamard_G_except!(V, Gs, n) + M_mttkrp = mttkrp!( mttkrp_bufs[n], A, U, @@ -306,10 +320,10 @@ function fit_cp_als( kr_work = mttkrp_kr_work2[n], ) if update_policy != :ls - _clamp_nonnegative!(G) + _clamp_nonnegative!(M_mttkrp) _nncp_mode_update!( U[n], - G, + M_mttkrp, V, denom_work[n]; nn_update = update_policy, @@ -317,28 +331,22 @@ function fit_cp_als( nnls_row_tol = nnls_row_tol, ) mul!(denom_work[n], U[n], V) - pg_sq += _projected_grad_sq_nonnegative(U[n], G, denom_work[n]) + pg_sq += _projected_grad_sq_nonnegative(U[n], M_mttkrp, denom_work[n]) u_sq += sum(abs2, U[n]) else - _als_mode_update!(U[n], G, V, transposed_work[n]) + _als_mode_update!(U[n], M_mttkrp, V, transposed_work[n]) end - _update_gram!(grams[n], U[n]) + _update_G!(Gs[n], U[n]) end - _cp_als_normalize_weights_and_grams!( - λ, - U, - grams, - normalization_policy, - update_policy, - ) + _cp_als_normalize_weights_and_Gs!(λ, U, Gs, normalization_policy, update_policy) _, _, rel_error = _cp_als_stats( A, normA2, λ, U, - grams; + Gs; mttkrp_method, mttkrp_buf = mttkrp_bufs[1], mttkrp_work = mttkrp_tmp_work[1], @@ -391,7 +399,7 @@ function fit_cp_als( pg_norm = _projected_grad_norm_nonnegative!( A, U, - grams, + Gs, V, denom_work, mttkrp_bufs, @@ -407,7 +415,7 @@ function fit_cp_als( normA2, λ, U, - grams; + Gs; mttkrp_method, mttkrp_buf = mttkrp_bufs[1], mttkrp_work = mttkrp_tmp_work[1], diff --git a/src/solvers/nncp_updates.jl b/src/solvers/nncp_updates.jl index b760bf1..b497ed4 100644 --- a/src/solvers/nncp_updates.jl +++ b/src/solvers/nncp_updates.jl @@ -47,17 +47,17 @@ end function _nncp_mu_mode_update!( U::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, denom::AbstractMatrix{T}, ) where {T<:AbstractFloat} floor = sqrt(eps(T)) - U .= max.(U .* G ./ (denom .+ floor), floor) + U .= max.(U .* M_mttkrp ./ (denom .+ floor), floor) return U end function _nncp_hals_mode_update!( U::StridedMatrix{T}, - G::StridedMatrix{T}, + M_mttkrp::StridedMatrix{T}, V::StridedMatrix{T}, work::StridedMatrix{T}, ) where {T<:AbstractFloat} @@ -67,7 +67,7 @@ function _nncp_hals_mode_update!( vkk = max(V[k, k], floor) for i in axes(U, 1) old = U[i, k] - new = max((G[i, k] - work[i, k] + old * vkk) / vkk, floor) + new = max((M_mttkrp[i, k] - work[i, k] + old * vkk) / vkk, floor) Δ = new - old U[i, k] = new if !iszero(Δ) @@ -82,7 +82,7 @@ end function _nncp_hals_mode_update!( U::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, V::AbstractMatrix{T}, work::AbstractMatrix{T}, ) where {T<:AbstractFloat} @@ -90,7 +90,7 @@ function _nncp_hals_mode_update!( d = max.(diag(V), floor) d_row = reshape(d, 1, :) mul!(work, U, V) - U .= max.((G .- work .+ U .* d_row) ./ d_row, floor) + U .= max.((M_mttkrp .- work .+ U .* d_row) ./ d_row, floor) return U end @@ -147,7 +147,7 @@ end function _nncp_nnls_mode_update!( U::StridedMatrix{T}, - G::StridedMatrix{T}, + M_mttkrp::StridedMatrix{T}, V::StridedMatrix{T}, work::StridedMatrix{T}; max_cd_sweeps::Int = 10, @@ -156,7 +156,7 @@ function _nncp_nnls_mode_update!( @inbounds for i in axes(U, 1) _nncp_nnls_row_update!( view(U, i, :), - view(G, i, :), + view(M_mttkrp, i, :), V, view(work, i, :); max_cd_sweeps, @@ -168,7 +168,7 @@ end function _nncp_nnls_mode_update!( U::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, V::AbstractMatrix{T}, work::AbstractMatrix{T}; max_cd_sweeps::Int = 10, @@ -179,7 +179,7 @@ function _nncp_nnls_mode_update!( d_row = reshape(d, 1, :) @inbounds for _ = 1:max_cd_sweeps mul!(work, U, V) - U_new = max.(U .- (work .- G) ./ d_row, floor) + U_new = max.(U .- (work .- M_mttkrp) ./ d_row, floor) max_delta = maximum(abs.(U_new .- U)) copyto!(U, U_new) max_delta <= row_tol * max(maximum(U), one(T)) && break @@ -189,7 +189,7 @@ end function _nncp_mode_update!( U::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, V::AbstractMatrix{T}, denom::AbstractMatrix{T}; nn_update, @@ -198,7 +198,7 @@ function _nncp_mode_update!( ) where {T<:AbstractFloat} return _nncp_mode_update!( U, - G, + M_mttkrp, V, denom, nn_update_policy(nn_update); @@ -209,7 +209,7 @@ end @inline function _nncp_mode_update!( U::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, V::AbstractMatrix{T}, denom::AbstractMatrix{T}, ::MultiplicativeNNUpdate; @@ -217,25 +217,25 @@ end nnls_row_tol::T = sqrt(eps(T)), ) where {T<:AbstractFloat} mul!(denom, U, V) - return _nncp_mu_mode_update!(U, G, denom) + return _nncp_mu_mode_update!(U, M_mttkrp, denom) end @inline function _nncp_mode_update!( U::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, V::AbstractMatrix{T}, denom::AbstractMatrix{T}, ::HALSNNUpdate; nnls_max_cd_sweeps::Int = 10, nnls_row_tol::T = sqrt(eps(T)), ) where {T<:AbstractFloat} - return _nncp_hals_mode_update!(U, G, V, denom) + return _nncp_hals_mode_update!(U, M_mttkrp, V, denom) end @inline function _nncp_mode_update!( U::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, V::AbstractMatrix{T}, denom::AbstractMatrix{T}, ::NNLSUpdate; @@ -244,7 +244,7 @@ end ) where {T<:AbstractFloat} return _nncp_nnls_mode_update!( U, - G, + M_mttkrp, V, denom; max_cd_sweeps = nnls_max_cd_sweeps, @@ -254,12 +254,12 @@ end @inline function _projected_grad_sq_nonnegative( U::AbstractMatrix{T}, - G::AbstractMatrix{T}, + M_mttkrp::AbstractMatrix{T}, denom::AbstractMatrix{T}, ) where {T<:AbstractFloat} floor = sqrt(eps(T)) active_floor = 10 * floor - grad = denom .- G + grad = denom .- M_mttkrp pg = ifelse.(U .<= active_floor, min.(grad, zero(T)), grad) return sum(abs2, pg) end @@ -271,16 +271,16 @@ function _projected_grad_norm_nonnegative!( V::AbstractMatrix{T}, denom_work::AbstractVector{<:AbstractMatrix{T}}, mttkrp_bufs::AbstractVector{<:AbstractMatrix{T}}, - mttkrp_tmp_work::AbstractVector{<:AbstractMatrix{T}}, - mttkrp_kr_work::AbstractVector{<:AbstractMatrix{T}}, - mttkrp_kr_work2::AbstractVector{<:AbstractMatrix{T}}; + mttkrp_tmp_work::AbstractVector, + mttkrp_kr_work::AbstractVector, + mttkrp_kr_work2::AbstractVector; mttkrp_method::Symbol = :auto, ) where {T<:AbstractFloat,N} sq = zero(T) u_sq = zero(T) for n = 1:N - _hadamard_gram_except!(V, grams, n) - G = mttkrp!( + _hadamard_G_except!(V, grams, n) + M_mttkrp = mttkrp!( mttkrp_bufs[n], A, U, @@ -290,9 +290,9 @@ function _projected_grad_norm_nonnegative!( kr_buf = mttkrp_kr_work[n], kr_work = mttkrp_kr_work2[n], ) - _clamp_nonnegative!(G) + _clamp_nonnegative!(M_mttkrp) mul!(denom_work[n], U[n], V) - sq += _projected_grad_sq_nonnegative(U[n], G, denom_work[n]) + sq += _projected_grad_sq_nonnegative(U[n], M_mttkrp, denom_work[n]) u_sq += sum(abs2, U[n]) end return sqrt(sq / max(u_sq, one(T))) diff --git a/src/tucker/sthosvd.jl b/src/tucker/sthosvd.jl index 8859e55..93a0609 100644 --- a/src/tucker/sthosvd.jl +++ b/src/tucker/sthosvd.jl @@ -18,9 +18,12 @@ by sequentially computing truncated SVDs mode by mode, projecting (shrinking) th # TuckerResult struct and show live in core/types.jl. Uses unfold_mode, mode_n_product from core/tensor_ops. """ - reconstruct(td::TuckerResult) reconstructs the tensor from Tucker decomposition + reconstruct(td::TuckerResult) - A = S ×₁ U₁ ×₂ U₂ ⋯ ×_d U_d +Reconstruct the dense tensor represented by a Tucker decomposition result. + +For a Tucker result with core `S` and factors `U_1, ..., U_d`, this returns +`S ×_1 U_1 ×_2 U_2 ... ×_d U_d`. """ function reconstruct(td::TuckerResult{T,N}) where {T,N} A = td.core diff --git a/test/basic_tests.jl b/test/basic_tests.jl index 0445ff7..c828ba4 100644 --- a/test/basic_tests.jl +++ b/test/basic_tests.jl @@ -2139,13 +2139,73 @@ end ) @test err_tuck_model <= err_diag_model + 1e-8 - # mttkrp: contract path should match explicit Khatri-Rao path + # mttkrp: direct path should match explicit Khatri-Rao path U3 = [randn(8, 3), randn(6, 3), randn(5, 3)] for mode = 1:3 G_kr = mttkrp(A, U3, mode; method = :khatri_rao) - G_ct = mttkrp(A, U3, mode; method = :contract) - @test G_ct ≈ G_kr atol = 1e-10 + G_dir = mttkrp(A, U3, mode; method = :direct) + @test G_dir ≈ G_kr atol = 1e-10 end + + A4 = randn(7, 5, 4, 3) + U4 = [randn(7, 2), randn(5, 2), randn(4, 2), randn(3, 2)] + for mode = 1:4 + G_kr = mttkrp(A4, U4, mode; method = :khatri_rao) + G_dir = mttkrp(A4, U4, mode; method = :direct) + @test G_dir ≈ G_kr atol = 1e-10 + end + + A5 = randn(4, 3, 2, 3, 2) + U5 = [randn(4, 2), randn(3, 2), randn(2, 2), randn(3, 2), randn(2, 2)] + for mode = 1:5 + G_kr = mttkrp(A5, U5, mode; method = :khatri_rao) + G_dir = mttkrp(A5, U5, mode; method = :direct) + @test G_dir ≈ G_kr atol = 1e-10 + end + + ws_direct3 = TensorKitchen.CPALSWorkspace(A, size(A), 3; mttkrp_method = :direct) + @test all(isnothing, ws_direct3.mttkrp_kr_work) + @test all(isnothing, ws_direct3.mttkrp_kr_work2) + @test all(x -> !isnothing(x), ws_direct3.mttkrp_tmp_work) + + ws_kr = TensorKitchen.CPALSWorkspace(A, size(A), 3; mttkrp_method = :khatri_rao) + @test all(x -> !isnothing(x), ws_kr.mttkrp_kr_work) + @test all(x -> !isnothing(x), ws_kr.mttkrp_kr_work2) + + ws_direct5 = TensorKitchen.CPALSWorkspace(A5, size(A5), 2; mttkrp_method = :direct) + @test all(isnothing, ws_direct5.mttkrp_kr_work) + @test all(isnothing, ws_direct5.mttkrp_kr_work2) + @test all(isnothing, ws_direct5.mttkrp_tmp_work) + + out_buf = zeros(size(A, 1), size(U3[1], 2)) + @test_throws ArgumentError TensorKitchen.mttkrp!( + out_buf, + A, + Matrix{Float64}[], + 1; + method = :direct, + ) + @test_throws DimensionMismatch TensorKitchen.mttkrp!( + out_buf, + A, + [U3[1], U3[2]], + 1; + method = :direct, + ) + @test_throws DimensionMismatch TensorKitchen.mttkrp!( + out_buf, + A, + [randn(7, 3), U3[2], U3[3]], + 1; + method = :direct, + ) + @test_throws DimensionMismatch TensorKitchen.mttkrp!( + out_buf, + A, + [U3[1], randn(6, 2), U3[3]], + 1; + method = :direct, + ) end @testset "utils: cross_component, build_cross_matrix, grad_lambda_cp, cp_rankr_cost_value, cross_term_gradU, gradU_column_cp" begin