Refactor join components and CP parameterizations

.github/workflows/format_check.yml

@@ -18,8 +18,7 @@ jobs:
       - uses: actions/checkout@v6
       - name: Install JuliaFormatter and format
         run: |
-          julia --project=. -e 'using Pkg; Pkg.instantiate()'
-          julia --project=. -e 'using JuliaFormatter; format(["./src", "./test"], verbose=true)'
+          julia -e 'using Pkg; Pkg.activate(temp=true); Pkg.add(PackageSpec(name="JuliaFormatter", version="2.10.1")); using JuliaFormatter; format(["./src", "./test"], verbose=true)'
       - name: Format check
         run: |
           julia -e '

Project.toml

@@ -14,10 +14,10 @@ RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2"
 
 [compat]
-JuliaFormatter = "2.8.5"
+JuliaFormatter = "=2.10.1"
 Manifolds = "0.11.28"
 ManifoldsBase = "2.3.5"
-Manopt = "0.5.37"
+Manopt = "0.6"
 ProgressMeter = "1.11.0"
 RecursiveArrayTools = "4.3"
 TensorOperations = "5.6"

src/api/btd.jl

@@ -63,6 +63,15 @@ _btd_uses_warm_start(::AbstractSolver, ::BTDALSWarmStartInit) = true
 _btd_should_polish(::ALSSolver, ::Integer) = false
 _btd_should_polish(::AbstractSolver, polish_n::Integer) = polish_n > 0
 
+function _reject_unsupported_btd_solver(solver_obj)
+    solver_obj isa LMSolver || return nothing
+    throw(
+        ArgumentError(
+            "BTD currently does not support LM refinement because the required Manopt operator path is not yet available for nested Tucker layouts. Use :rgd, :rcg, :lbfgs, :als, or :btd_tsd instead.",
+        ),
+    )
+end
+
 function _btd_warm_start_result(
     model::JoinModel{T,<:BTDBackend},
     backend::BTDBackend,
@@ -159,7 +168,10 @@ refines it. Returns a [`BTDResult`](@ref).
     - `:als`: Alternating least squares.
     - `:rcg`: Riemannian conjugate gradient.
     - `:lbfgs`: Limited-memory quasi-Newton refinement.
-    - `:lm`: Levenberg-Marquardt refinement.
+    - `:btd_tsd`: Blockwise tangent-subspace descent for BTD.
+
+`solver = :lm` is currently not supported for BTD because the required Manopt
+LM operator path is not yet available for nested Tucker layouts.
 
 ## Extended Options
 
@@ -234,6 +246,7 @@ function btd(
     kwargs...,
 ) where {T<:AbstractFloat,N}
     solver_obj = _solver_object(solver, stepsize; kwargs...)
+    _reject_unsupported_btd_solver(solver_obj)
     solver_sym = _btd_solver_symbol(solver_obj)
     init_resolved = _resolve_btd_init(init, solver_obj)
     init_eff =

src/btd/model.jl

@@ -31,8 +31,18 @@ Backend state for block-term decomposition as a sum of Tucker blocks. Stores
 the target tensor, product manifold, reusable work buffers, and per-block
 ambient reconstruction buffers used by cost, gradient, and ALS routines.
 """
-struct BTDBackend{T,N,MT<:Tuple,A<:AbstractArray{T,N},V,MP<:ProductManifold,I,C} <:
-       AbstractJoinBackend
+struct BTDBackend{
+    T,
+    N,
+    CT<:Tuple,
+    MT<:Tuple,
+    A<:AbstractArray{T,N},
+    V,
+    MP<:ProductManifold,
+    I,
+    C,
+} <: AbstractJoinBackend
+    components::CT
     manifolds::MT
     r::Int
     # Preserve the target array/backend so BTD shares the generic join storage behavior.
@@ -60,7 +70,7 @@ model_exact_join_basis_function(model::JoinModel{<:AbstractFloat,<:BTDBackend})
     (M, p) -> begin
         backend = model.backend
         residual = _join_residual!(backend, p)
-        _join_basis_project(backend.manifolds, p, residual)
+        _join_basis_project(backend.components, p, residual)
     end
 
 
@@ -71,22 +81,19 @@ function _btd_sequential_tucker_init(model::JoinModel{<:AbstractFloat,<:BTDBacke
     residual = copy(backend.target)
     parts = Vector{Manifolds.TuckerPoint{eltype(backend.target)}}(undef, backend.r)
     for k = 1:backend.r
-        pk = _manifold_init(backend.manifolds[k], residual, init)
+        component = _backend_component(backend, k)
+        Mk = _backend_manifold(backend, k)
+        pk = _component_init(component, residual, init)
         parts[k] = pk
-        _subtract_ambient_tensor!(
-            residual,
-            backend.manifolds[k],
-            pk,
-            backend.component_bufs[k],
-        )
+        _subtract_ambient_tensor!(residual, Mk, pk, backend.component_bufs[k])
     end
     return ArrayPartition(parts...)
 end
 
 function _btd_block_ranks_by_mode(backend::BTDBackend{T,N}) where {T,N}
     ranks = Vector{NTuple{N,Int}}(undef, backend.r)
     for b = 1:backend.r
-        M = backend.manifolds[b]
+        M = _backend_manifold(backend, b)
         M isa Manifolds.Tucker || throw(
             ArgumentError(
                 "BTD HOSVD multistart expects Tucker manifolds, got $(typeof(M)) at block $b.",
@@ -160,7 +167,7 @@ function _btd_hosvd_split_candidate(
         parts[b] = pk
         _subtract_ambient_tensor!(
             residual,
-            backend.manifolds[b],
+            _backend_manifold(backend, b),
             pk,
             backend.component_bufs[b],
         )
@@ -190,7 +197,7 @@ function initial_point(
     init_sym = _builtin_initializer_symbol(init)
     if init_sym == :random
         parts = ntuple(
-            k -> _manifold_init(backend.manifolds[k], backend.target, init),
+            k -> _component_init(_backend_component(backend, k), backend.target, init),
             backend.r,
         )
         return ArrayPartition(parts...)
@@ -282,6 +289,9 @@ function rgrad(model::JoinModel{<:AbstractFloat,<:BTDBackend}, p)
     parts = point_parts(p)
     _check_parts_len(parts, backend.r, "BTD rgrad")
     eg = point_parts(_btd_egrad(backend, p))
-    vals = ntuple(k -> egrad_to_rgrad(backend.manifolds[k], parts[k], eg[k]), backend.r)
+    vals = ntuple(
+        k -> egrad_to_rgrad(_backend_manifold(backend, k), parts[k], eg[k]),
+        backend.r,
+    )
     return wrap_like_point(p, vals)
 end

src/core/model.jl

@@ -1,6 +1,17 @@
 # core/model.jl — Top-level decomposition model interface
 export AbstractDecompositionModel,
-    manifold, initial_point, egrad, rgrad, supports_rgrad, tensor, cost, post_step!
+    manifold,
+    initial_point,
+    egrad,
+    rgrad,
+    supports_rgrad,
+    tensor,
+    cost,
+    post_step!,
+    residual,
+    differential_action,
+    differential_action!,
+    adjoint_action
 """
     AbstractDecompositionModel{T}
 
@@ -113,3 +124,49 @@ end
 function tensor(model::AbstractDecompositionModel)
     error("tensor not implemented for $(typeof(model))")
 end
+
+function residual(model::AbstractDecompositionModel, p)
+    error("residual not implemented for $(typeof(model))")
+end
+
+function differential_action!(out::AbstractVector, model::AbstractDecompositionModel, p, X)
+    error("differential_action! not implemented for $(typeof(model))")
+end
+
+function differential_action(model::AbstractDecompositionModel, p, X)
+    T = eltype(tensor(model))
+    out = Vector{T}(undef, length(tensor(model)))
+    differential_action!(out, model, p, X)
+    return out
+end
+
+function adjoint_action(
+    model::AbstractDecompositionModel,
+    p,
+    a::AbstractVector;
+    basis = ManifoldsBase.DefaultOrthonormalBasis(),
+)
+    M = manifold(model)
+    d = manifold_dimension(M)
+    length(a) == length(tensor(model)) || throw(
+        DimensionMismatch(
+            "adjoint_action expected ambient vector of length $(length(tensor(model))) for $(typeof(model)), got $(length(a)).",
+        ),
+    )
+    T = _scalar_eltype(p)
+    coeff = zeros(T, d)
+    e_j = zeros(T, d)
+    col = Vector{T}(undef, length(a))
+    @inbounds for j = 1:d
+        fill!(e_j, zero(T))
+        e_j[j] = one(T)
+        Xj = ManifoldsBase.get_vector(M, p, e_j, basis)
+        differential_action!(col, model, p, Xj)
+        coeff[j] = dot(col, a)
+    end
+    return ManifoldsBase.get_vector(M, p, coeff, basis)
+end
+
+function adjoint_action(model::AbstractDecompositionModel, p, a::AbstractArray; kwargs...)
+    return adjoint_action(model, p, vec(a); kwargs...)
+end