Update docs

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
 
-
 <br>
 
+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`.

docs/make.jl

@@ -1,7 +1,7 @@
 using Documenter, TensorKitchen
 
 makedocs(
-    sitename = "Documentation",
+    sitename = "TensorKitchen.jl",
     warnonly = true,
     pages = [
         "Home" => "index.md",

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`

docs/src/index.md

@@ -1,5 +1,6 @@
-# TensorKitchen.jl Documentation
+# TensorKitchen.jl: tensor decompositions in Julia
 
+<img src="logo_text_transparent.png" width="450px">
 
 **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/).
 

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.

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))

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

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."))