Add Riemannian Levenberg-Marquardt solver via Manopt.

src/api/approx.jl

@@ -127,11 +127,12 @@ For the generic join path:
     - `rgd_fixed`: Riemannian gradient descent with fixed step size
     - `rcg`: Riemannian conjugate gradient
     - `lbfgs`: Limited-memory quasi-Newton
+    - `lm`: Levenberg-Marquardt on residual/Jacobian least squares
 
 ##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 random initial point and then use manifold solvers for refinement.
-* For generic mixed joins, use manifold solvers such as `:rgd`, `:rcg`, or `:lbfgs`.
+* For generic mixed joins, use manifold solvers such as `:rgd`, `:rcg`, `:lbfgs`, or `:lm`.
 """
 function _approx_manifold_collection(
     dispatch::AutoApproxDispatch,

src/api/btd.jl

@@ -159,6 +159,7 @@ refines it. Returns a [`BTDResult`](@ref).
     - `:als`: Alternating least squares.
     - `:rcg`: Riemannian conjugate gradient.
     - `:lbfgs`: Limited-memory quasi-Newton refinement.
+    - `:lm`: Levenberg-Marquardt refinement.
 
 ## Extended Options
 

src/api/cpd.jl

@@ -134,7 +134,7 @@ function _component_energy_summary(energies::AbstractVector{<:Real})
         top1 = ordered[1],
         top2 = sum(@view ordered[1:min(2, r)]),
         top3 = sum(@view ordered[1:min(3, r)]),
-        effective = hhi > 0 ? inv(hhi) : NaN,
+        effective = hhi > 0 ? one(hhi) / hhi : NaN,
         argmax_component = argmax(shares),
     )
 end
@@ -579,13 +579,13 @@ end
 function _validate_cpd_solver_supported(solver::AbstractSolver)
     throw(
         ArgumentError(
-            "Unsupported CPD solver $(typeof(solver)). Use :als, :rgd, :rgd_fixed, :rcg, or :lbfgs.",
+            "Unsupported CPD solver $(typeof(solver)). Use :als, :rgd, :rgd_fixed, :rcg, :lbfgs, or :lm.",
         ),
     )
 end
 
 _validate_cpd_solver_supported(
-    ::Union{ALSSolver,RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver},
+    ::Union{ALSSolver,RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver,LMSolver},
 ) = nothing
 
 function _validate_cpd_solver_options(
@@ -750,7 +750,7 @@ end
 
 function _cpd_manifold_grad_tol(
     model::JoinModel{<:AbstractFloat,<:CPDBackend},
-    solver::Union{RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver},
+    solver::Union{RGDSolver,RGDFixedSolver,RCGSolver,LBFGSSolver,LMSolver},
     tol::Real,
 )
     return tol
@@ -992,6 +992,7 @@ If `r` is omitted, uses the smallest tensor mode as a heuristic rank.
     - `rgd_fixed`: Riemannian gradient descent with fixed step size
     - `rcg`: Riemannian conjugate gradient
     - `lbfgs`: Limited-memory Riemannian quasi-Newton
+    - `lm`: Levenberg-Marquardt using residual/Jacobian least squares
     - `als`: Alternating Least Squares
 
 ## Extended Options

src/backend.jl

@@ -8,6 +8,7 @@ include("solvers/rgd.jl")
 include("solvers/btd_tsd.jl")
 include("solvers/rcg.jl")
 include("solvers/lbfgs.jl")
+include("solvers/lm.jl")
 include("results/reconstruct.jl")
 include("results/rel_error.jl")
 include("solvers/solve_dispatch.jl")

src/core/model.jl

@@ -1,5 +1,6 @@
 # core/model.jl — Top-level decomposition model interface
-export AbstractDecompositionModel, rgrad, supports_rgrad, tensor, cost, post_step!
+export AbstractDecompositionModel,
+    manifold, initial_point, egrad, rgrad, supports_rgrad, tensor, cost, post_step!
 """
     AbstractDecompositionModel{T}
 

src/core/types.jl

@@ -446,7 +446,7 @@ function LinearAlgebra.cond(R::CPDResult{T}) where {T<:AbstractFloat}
     U = hcat(Us...)
     s = svdvals(U)
     smin = minimum(s)
-    return iszero(smin) ? T(Inf) : inv(smin)
+    return iszero(smin) ? T(Inf) : one(T) / smin
 end
 
 """

src/core/unpack_points.jl

@@ -1,6 +1,11 @@
 # core/unpack_points.jl — Point unpacking and legacy/vector interop
 export pack_point_rank1,
-    unpack_point_rank1, pack_point_rankr, unpack_point_rankr, unpack_point_rankr_components
+    unpack_point_rank1,
+    pack_point_rankr,
+    unpack_point_rankr,
+    unpack_point_rankr_components,
+    canonical_to_joinpoint,
+    joinpoint_to_canonical
 
 function unpack_rankr_native(p, dims::NTuple{N,Int}, r::Int) where {N}
     parts = normalize_rankr_native_point(p, dims, r)
@@ -72,6 +77,43 @@ function unpack_rankr_join(p, dims::NTuple{N,Int}, r::Int) where {N}
     return λ, U
 end
 
+"""
+    canonical_to_joinpoint(λ, U)
+    canonical_to_joinpoint(p_canonical, dims, r)
+
+Convert a CPD point from canonical factor-matrix storage to the native
+rank-`r` Segre join point layout used by generic `JoinModel((Segre, ...), A)`.
+
+The conversion preserves the represented tensor but may renormalize component
+gauges the same way `pack_rankr_native` does.
+"""
+function canonical_to_joinpoint(
+    λ::AbstractVector{T},
+    U::Vector{<:AbstractMatrix{T}},
+) where {T<:AbstractFloat}
+    r = length(λ)
+    return pack_rankr_native(λ, U, r)
+end
+
+function canonical_to_joinpoint(p, dims::NTuple{N,Int}, r::Int) where {N}
+    λ, U = unpack_rankr_canonical(p, dims, r)
+    return canonical_to_joinpoint(λ, U)
+end
+
+"""
+    joinpoint_to_canonical(p_join, dims, r)
+
+Convert a native Segre join point layout back to the canonical CPD point
+layout `(λ, (u₁¹, …, uᵣ¹), …, (u₁ᴺ, …, uᵣᴺ))`.
+
+The conversion preserves the represented tensor but may renormalize component
+gauges the same way `pack_rankr_canonical` does.
+"""
+function joinpoint_to_canonical(p, dims::NTuple{N,Int}, r::Int) where {N}
+    λ, U = unpack_rankr_native(p, dims, r)
+    return pack_rankr_canonical(λ, U, r)
+end
+
 function pack_point_rank1(λ::T, U::Vector{Vector{T}}) where {T<:AbstractFloat}
     parts = Vector{Vector{T}}(undef, length(U) + 1)
     parts[1] = T[λ]

src/cpd/core/cp_cost.jl

@@ -69,7 +69,7 @@ end
 end
 
 @inline function _softplus_derivative(x::Real)
-    x >= 0 ? inv(one(x) + exp(-x)) : begin
+    x >= 0 ? one(x) / (one(x) + exp(-x)) : begin
         ex = exp(x)
         ex / (one(x) + ex)
     end