Smooth Reparameterizations of Functions on Simplicial Product Spaces: Applications to Probabilistic Tensor Decomposition and Functional Data Registration

JOTA manuscript No. (will be inserted by the editor)

Smooth Reparameterizations of Functions on Simplicial Product Spaces: Applications to Probabilistic Tensor Decomposition and Functional Data Registration Shashwat Kumar · Arafat Rahman · Anuj Srivastava · P.-A. Absil

arXiv:2608.02576v1 [cs.LG] 3 Aug 2026

Received: date / Accepted: date

Abstract We consider optimization problems defined on product spaces of simplices. Examples of this class of problems include learning low-rank discrete multivariate probability distributions via simplex constrained tensor decomposition and performing functional data registration under the Square Root Velocity Function (SRVF) representation. In this work, we demonstrate the feasibility of replacing the product simplex with a smooth, elementwise strictly convex reparameterization, resulting in an unconstrained optimization problem on a manifold. We show that performing such a reparameterization results in the second order Karush–Kuhn–Tucker (KKT) points on the smooth manifold being mapped to the weak second order KKT points on the product simplex. This leads to a Riemannian Gradient Descent (RGD) algorithm for solving the reparameterized problem, which outperforms Projected Gradient Descent (PGD), and provides a more faithful representation of the original function shapes while performing curve registration. Keywords Riemannian Optimization · Tensor Decomposition · Functional Data Registration · Hadamard Parametrization Mathematics Subject Classification (2020) 65K05 · 53B21 · 70G45 Shashwat Kumar Department of Biomedical Engineering, Johns Hopkins University, USA E-mail: skuma118@jh.edu Arafat Rahman Department of Systems and Information Engineering, University of Virginia, USA E-mail: jgh6ds@virginia.edu Anuj Srivastava Department of Applied Mathematics and Statistics, Johns Hopkins University, USA E-mail: anuj.srivastava@jhu.edu P.-A. Absil ICTEAM Institute, UCLouvain, 1348 Louvain-la-Neuve, Belgium E-mail: pa.absil@uclouvain.be

2

Shashwat Kumar et al.

1 Introduction Several optimization problems involve a twice differentiable function f being optimized over a feasible set that is a product space of simplices and Euclidean components: min

v∈Fsimplex

f (v),

Fsimplex = ∆n1 × · · · × ∆nL × RnL+1 × · · · × RnL+Q .

(Psimplex )

Here, ∆n = {x ∈ Rn | x ≥ 0, 1T x = 1} denotes the standard (n − 1)simplex embedded in Rn . Such problems often naturally arise in diverse domains, such as (i) simplex-constrained matrix factorization, where the elements of the simplices may for example represent unknown mixture proportions of pure minerals to be estimated from hyperspectral images [1, 18, 22, 30]; (ii) probabilistic tensor decomposition, where simplices are used to learn probability distributions [17], perform feature selection [4], or minimize the smallest element of a tensor [26]; and (iii) shape analysis, where the Square Root Velocity Function (SRVF) representation leverages simplicial parameterizations to temporally align functional data [27, 28]. A common technique for solving optimization problems with such constraints involves strategies like Projected Gradient Descent (PGD) which involve moving along the direction of the gradient followed by projection on Fsimplex [18]. As an alternative for solving (Psimplex ), several papers have focused on solving by Riemannian optimization algorithms a smooth reparameterized version of the problem: min

θ∈Fsmooth

f (φ(θ)),

(Psmooth )

where Fsmooth is a manifold and φ : Fsmooth → Fsimplex is surjective. For instance, the Hadamard reparameterization φ : Sn−1 → ∆n , defined by φ(θ) = θ ⊙ θ, where Sn−1 := {θ ∈ Rn : ∥θ∥2 = 1} denotes the unit sphere in Rn , has been used in linear programming [12], sparse recovery [29], and highdimensional linear regression [32]. Products of Hadamard maps φ : Sn−1 ×

n−1 n n ··· × S → ∆ × · · · × ∆ , φ(θ) = θ1 ⊙ θ1 , . . . , θL ⊙ θL have been used to solve Simplex Constrained Matrix Factorization in [11] and [15]. Other works also add Euclidean components for use in shape registration [28] φ : Sn−1 × ·

· · × Sn−1 × Rn → ∆n × · · · × ∆n × Rn , φ(θ) = θ1 ⊙ θ1 , . . . , θL ⊙ θL , θL+1 . Non-Hadamard parameterizations φ : Rn × · · · × Rn → ∆n × · · · × θ1 θL

∆n , φ(θ) = 1Te eθ1 , . . . , 1Te eθL have also been used in probabilistic tensor factorizations [26]. These algorithms can solve (Psmooth ) by moving along the manifold Fsmooth without requiring projection back to Fsimplex . A recent body of work has studied the critical points of (Psimplex ) and (Psmooth ) under the Hadamard map, φ : Sn−1 → ∆n , φ(θ) = θ ⊙ θ, showing that second order KKT points for (Psmooth ) imply first order [20] and (weak) second-order criticality [21] for (Psimplex ), respectively. However, several of the

Title Suppressed Due to Excessive Length

3

applications we described above involve more general feasible sets and transformations, and the existing results do not directly extend to these settings. Furthermore, while the Hadamard map has recieved a lot of attention, it is not known if there are other more general classes of transformations which preserve KKT points. Contributions: – We generalize the results of Li et al. [21] from Hadamard parameterization to a more general class of smooth transformations ψ(θ), with the following properties: ψ : R → R≥0 , ψ(0) = 0,

ψ ′′ > 0.

(1)

Apart from the Hadamard map, this includes polynomials with even degree, positive coefficients, and no constant term, e.g., ψ(θ) = aθ4 + bθ2 , a, b > 0, and other functions such as ψ(θ) = cosh θ − 1. – We also generalize their results from Fsimplex = ∆n to Fsimplex = ∆n1 × · · · × ∆nL × RnL+1 × · · · × RnL+Q , allowing this result to be applied to problems like matrix factorization [11], tensor factorization [26], and functional data registration [28]. This makes it possible to tackle the problems by Riemannian optimization algorithms [2, 9]. – We distinguish a weak and a strong version of the second-order KKT conditions for (Psimplex ). We show that in Section 4 that the weak—but not the strong—second-order KKT points of (Psimplex ) correspond to the (weak, equivalently strong) second-order KKT points of the parameterized formulation. – We provide numerical results on two problems in particular: a) simplex constrained tensor factorization and b) shape registration. We introduce a Riemannian Gradient Descent (RGD) algorithm which outperforms PGD for simplex constrained tensor factorization problems and achieves smoother and more accurate registration for the shape registration problem.

2 Reparameterization of (Psimplex ) We reparameterize (Psimplex ) under transformations of the form v = φ(θ) = [ψ(θ1 ), . . . , ψ(θL ), θL+1 , . . . , θL+Q ] with ψ as in (1). This leads to a reparameterized problem: min

θ∈Fsmooth

f (φ(θ)),

Fsmooth = {θ ∈ Rn1 × · · · × RnL+Q | 1T ψ(θl ) = 1,

l = 1, . . . , L}. (P×smooth )

Here, function ψ applied to a vector is meant elementwise: (ψ(θl ))i = ψ(θli ). As for ψ ′ and ψ ′′ , they are the elementwise first and second-order derivative,

4

Shashwat Kumar et al.

Fig. 1 a) Example of a function f (v) = v T Av on Fsimplex . b) Reparameterizing f (v) on a smooth manifold Fsmooth as f (φ(θ)), and c) its Legendre transform f (φ∗ (θ∗ )). Darker colors correspond to lower values of f . At the lower right corner of the simplex, all the equations in the KKT conditions (4) can be satisfied except βl ≥ 0, hence the lower right corner is not a first-order KKT point; its counterpart in Fsmooth is the point with the lightest color which, being a maximizer subject to a submanifold constraint, is a KKT point.

respectively. Fig. 1 shows an example of this type of reparameterization. Two important properties of ψ immediately stand out: Remark 2.1 ψ ′ (θ) = 0 if and only if θ = 0. Indeed, ψ is smooth and 0 is a minimizer, hence ψ ′ (0) = 0; and if θ ̸= 0, then the mean value theorem applied to ψ ′ yields ψ ′ (θ) = ψ ′′ (θ̃)θ, which is nonzero since ψ ′′ > 0. Remark 2.2 ψ(θ) = 0 if and only if θ = 0. Indeed, the condition is sufficient by definition of ψ; for the converse, observe that Taylor’s theorem with the remainder in Lagrange form yields ψ(θ) = 12 ψ ′′ (θ̃)θ2 > 0 for all θ ̸= 0. Remark 2.3 The range of ψ is R≥0 . Indeed, since ψ ′′ > 0, the function ψ is convex and nonconstant, and therefore it cannot be bounded from above, as every convex function on R that is bounded from above is constant [23, § 1.5, Exercise 2]. Since ψ(0) = 0, ψ ≥ 0, and ψ is continuous, the intermediate value theorem yields Im(ψ) = R≥0 . Proposition 2.1 The feasible set Fsmooth is a Riemannian manifold. Proof Let us show that Fsmooth is a Cartesian product of embedded submanifolds. In view of [2, Proposition 3.3.3], it suffices to show that the differential of c : Rnl → R : θl 7→ 1T ψ(θl ) is surjective, namely nonzero, at all θl that satisfy the simplex constraint c(θl ) = 1. The Jacobian of c is the row vector (ψ ′ (θl ))T . In view of the above remarks, the simplex constraint c(θl ) = 1 implies that at least one component of vector θl is nonzero, which ensures that (ψ ′ (θl ))T is nonzero. ⊔ ⊓ We now study relationships between critical points of (Psimplex ) and (P×smooth ).

Title Suppressed Due to Excessive Length

5

3 KKT Analysis Constructing the Lagrangians of (Psimplex ) and (P×smooth ), we get, L1 (v, λ, β) = f (v) +

L X

λl (1T vl − 1) − βlT vl ,

(2)

l=1

s

L2 (θ, λ ) = f (φ(θ)) +

L X

λs (1T ψ(θl ) − 1).

(3)

l=1

The first-order KKT conditions of (Psimplex ) take the form ∇vl f (v) + λl 1 − βl = 0,

l = 1, . . . , L,

(4a)

βl ⊙ vl = 0,

l = 1, . . . , L,

(4b)

vl ≥ 0, βl ≥ 0,

l = 1, . . . , L,

(4c)

l = 1, . . . , L,

(4d)

q = L + 1, . . . , L + Q,

(4e)

T

1 vl = 1, ∇vq f (v) = 0, and those of (P×smooth ) read

ψ ′ (θl ) ⊙ ∇vl f (v) + λs ψ ′ (θl ) = 0, 1T ψ(θl ) = 1, ∇θq f (φ(θ)) = 0,

l = 1, . . . , L,

(5a)

l = 1, . . . , L,

(5b)

q = L + 1, . . . , L + Q.

(5c)

They reduce to those in Li et al. [21, §4] when L = 1, Q = 0, and ψ(θ) = θ ⊙ θ. Let v ⋆ and θ⋆ be feasible points. Assume that (v ⋆ , λ⋆ , β ⋆ ) and (θ⋆ , λs,⋆ ) satisfy the first-order KKT conditions (4) and (5), respectively. The weak second-order KKT conditions for (Psimplex ) are given by   ∆v T ∇2v L1 (v ⋆ , λ⋆ , β ⋆ ) ∆v ≥ 0, ∆v = . . . , ∆vl , . . . , ∆vq , . . . , (6a) ∆vli = 0

1T ∆vl = 0,

l = 1, . . . , L,

(6b)

if (vl⋆ )i = 0,

l = 1, . . . , L,

(6c)

and the (weak, equivalently strong) second-order KKT conditions for (P×smooth ) are given by   ∆θT ∇2θ L2 (θ⋆ , λs,⋆ ) ∆θ ≥ 0, ∆θ = . . . , ∆θl , . . . , ∆θq , . . . , (7a) 1T (ψ ′ (θl⋆ ) ⊙ ∆θl ) = 0,

l = 1, . . . , L.

(7b)

Remark 3.1 Equations (5) express that the Riemannian gradient is zero for the Riemannian optimization problem (P×smooth ). Likewise, (7) express that the Riemannian Hessian is positive semidefinite.

6

Shashwat Kumar et al.

Remark 3.2 Condition (6) represents a weak form of the second-order necessary conditions. A weak second-order KKT point is a point v ⋆ for which there exist λ⋆ and β ⋆ such that (4) and (6) are satisfied. Note that this is distinct from the classical second-order necessary conditions for optimality. While Bertsekas et al. [6] and Gill et al. [13] present similar weak conditions, Bonnans et al. [8], Nocedal et al. [24], and Aragón et al. [5] present the strong conditions. For the strong version, we replace ∆vli = 0

if (vl⋆ )i = 0

with ∆vli = 0

if (vl⋆ )i = 0 and (βl⋆ )i > 0,

∆vli ≥ 0

if (vl⋆ )i = 0 and (βl⋆ )i = 0.

For (P×smooth ), the weak and strong conditions coincide because there are no inequality constraints. The weak/strong terminology follows the one in [14, Definition 2.1].

4 Main Results The Hessian from the second-order KKT of (Psimplex ) is given by: ∇2v L1 (v, λ, β) = ∇2v f (v).

(8)

Let θa and θb be two blocks of variables with dimensions na and nb respectively. Let g(θ) be a differentiable function. Define matrix ∇2θa θb g(θ) ∈ Rna ×nb by:
∂ 2 g(θ) . ∇2θa θb g(θ) ij = ∂θai ∂θbj Next, consider the Hessian of (P×smooth ). Here, similar to Li et al. [21, Equation 13], the Hessian is the sum of a diagonal term and a generally non-diagonal term involving the Hessian of f . However, our formulation includes multiple simplex and Euclidean blocks, leading to mixed second-derivative blocks, so the results of Li et al. [21] do not apply directly. For l, m ∈ {1, . . . , L} with l ̸= m and q ∈ {L + 1, . . . , L + Q}, we obtain ∇2θq θq L2 = ∇2vq vq f, ∇2θq θl L2 = ∇2vq vl f diag(ψ ′ (θl )), ∇2θl θq L2 = diag(ψ ′ (θl ))∇2vl vq f, ∇2θl θl L2 = diag(ψ ′′ (θl ) ⊙ ∇vl f + λsl ψ ′′ (θl )) + diag(ψ ′ (θl ))∇2vl vl f diag(ψ ′ (θl )), ∇2θm θl L2 = diag(ψ ′ (θm ))∇2vm vl f diag(ψ ′ (θl )). (9)

Title Suppressed Due to Excessive Length

7

Using (9), we can rewrite ∇2θ L2 as the sum of a diagonal matrix and a change of basis Hessian term   .. .    diag(ψ ′′ (θl ) ⊙ ∇vl f + λsl ψ ′′ (θl ))      .. ∇2θ L2 = +diag(φ′ (θ)) ∇2v f diag(φ′ (θ)). .     0q   .. . (10)   .. .    ψ ′ (θl )   
  . ′ .. Here, diag φ (θ) =   is the Jacobian of φ.     Iq   .. . We now proceed to show the following: Theorem 4.1 If v ⋆ is a first-order KKT point of (Psimplex ) then all θ⋆ satisfying v ⋆ = φ(θ⋆ ) are first-order KKT points of (P×smooth ). Proof We extend Li et al. [21, Equations 16–17] to the reparameterization proposed in (P×smooth ). Let us start with the KKT zero-gradient conditions with respect to each vector lying in Fsimplex (4a), ∇vl f + λl 1 − βl = 0,

l = 1, . . . , L.

(11)

Hadamarding with ψ ′ (θl ) and setting λl = λsl , ∇vl f ⊙ ψ ′ (θl ) + λsl ψ ′ (θl ) − βl ⊙ ψ ′ (θl ) = 0,

l = 1, . . . , L.

(12)

From the first-order KKT condition on the simplex (4b), we have βl ⊙ vl = 0,

l = 1, . . . , L.

(13)

This implies the following two cases: – If vli = 0, then θli = 0, leading to ψ ′ (θl )i = 0 (see Remarks 2.1, 2.2), and consequently βli ψ ′ (θl )i = 0. – If vli ̸= 0, then βli = 0, which also ensures βli ψ ′ (θl )i = 0. Thus, we conclude that βl ⊙ vl = 0 =⇒ βl ⊙ ψ ′ (θl ) = 0,

l = 1, . . . , L.

(14)

Plugging this into the equation above satisfies the zero-gradient conditions of (5a): ∇vl f ⊙ ψ ′ (θl ) + λsl ψ ′ (θl ) = 0, l = 1, . . . , L. (15)

8

Shashwat Kumar et al.

Next, the simplex constraint (5b) yields

1T vl = 1 =⇒ 1T ψ(θl ) = 1,

l = 1, . . . , L.

(16)

Finally, the Euclidean gradient conditions (5c) can be rewritten as

∇vq L1 = 0 =⇒ ∇vq f = 0 =⇒ ∇θq f = 0 =⇒ ∇θq L2 = 0.

(17)

Thus, a first-order KKT point on the product simplex implies a first-order KKT point on the smooth manifold. ⊓ ⊔

Theorem 4.2 If v ⋆ is a weak second-order KKT point of (Psimplex ) then all θ⋆ satisfying v ⋆ = φ(θ⋆ ) are second-order KKT points of (P×smooth ).

Proof We now use our Hessian formula in (10) to extend Li et al. [21, Equations 18-23] to show that second-order KKT points on Fsimplex implies second-order KKT points on Fsmooth . We start with an arbitrary vector in the tangent space of Fsmooth :

  ∆θ = . . . ∆θq . . . ∆θl . . . ,

∆θl s.t. 1T (ψ ′ (θl ) ⊙ ∆θl ) = 0,

l = 1, . . . , L. (18) Since ψ ′ (θl )⊙∆θl = diag(ψ ′ (θl ))∆θl = ∆vl , the corresponding vector ∆vl after applying the Jacobian to ∆θl satisfies 1T ∆vl = 0, l = 1, . . . , L. Furthermore, if vli = 0, then θli = 0 and ψ ′ (θl )i = 0 (Remark 2.1, 2.2) and thus ∆vli = 0.

Title Suppressed Due to Excessive Length

9

∆vl is thus a valid tangent vector on the corresponding simplex. ∆θT ∇2θ L2 ∆θ = 

 .. .    ∆θq      ..   . +     ∆θl    .. .. . . 

..

 .  0q    . .. . . . ∆θq . . . ∆θl . . .    diag(ψ ′′ (θl ) ⊙ ∇v f + λs ψ ′′ (θl )) l l 

∆θT diag(φ′ (θ))∇2v f diag(φ′ (θ))∆θ 



..

.   ′′ s = ∆θ1 . . . ∆θL   diag(ψ (θl ) ⊙ (∇vl f + λl 1))

..

 ∆θ1  .   . + .  .

∆θL

∆v T ∇2v f ∆v 

..

  . ′′ = ∆θ1 . . . ∆θL   diag(ψ (θl ) ⊙ βl )

=

nl L X X



 ∆θ1  .    .  + ∆v T ∇2 f ∆v .  .. ∆θ L .

ψ ′′ (θl )i βli (∆θli )2 + ∆v T ∇2 f ∆v.

(19)

l=1 i=1

Note that βli ≥ 0 (4c), ψ ′′ (θl )i > 0 (1) and ∆v T ∇2 f ∆v is nonnegative (6a), implying that (19) is nonnegative. ⊓ ⊔ It is interesting to note that Theorem 4.1 does not imply bidirectional criticality. We need an additional assumption of second-order criticality of (P×smooth ) in order to prove first-order criticality of (Psimplex ). Recall that the weak and strong second-order KKT points coincide for (P×smooth ). Theorem 4.3 If θ⋆ is a second-order KKT point of (P×smooth ) then v ⋆ = φ(θ⋆ ) is a first order KKT point of (Psimplex ). Proof Similar to Li et al. [21, Equations 24-31], we show that second-order KKT on smooth manifold implies first-order KKT point on Fsimplex . From KKT conditions on the Smooth manifold, we have: ψ ′ (θl ) ⊙ ∇vl f + λsl ψ ′ (θl ) = 0,

l = 1, . . . , L.

(20)

Since the corresponding point on Fsimplex is not a stationary point, the gradient conditions on the Simplex are: ∇vl f + λl 1 − βl = gl ,

l = 1, . . . , L.

(21)

10

Shashwat Kumar et al.

Absorbing gl in βl , these conditions can be rewritten as: ∇vl f + λl 1 − β̂l = 0,

(22)

l = 1, . . . , L.

Hadamarding equation (22) with ψ ′ (θl ) and setting λl = λsl ψ ′ (θl ) ⊙ ∇vl f + λsl ψ ′ (θl ) − ψ ′ (θl ) ⊙ β̂l = 0 =⇒ ψ ′ (θl ) ⊙ β̂l = 0,

(23) l = 1, . . . , L.

Furthermore because our choice of parameterization, ψ is 0 only at 0 (2.2). ψ ′ (θl ) ⊙ β̂l = 0 =⇒ θl ⊙ β̂l = 0 =⇒ vl ⊙ β̂l = 0,

(24)

l = 1, . . . , L.

Note that because of the nature of the square mapping, the positivity and summing up to 1, KKT conditions of (4) are already satisfied. As long as we can show that β̂l ≥ 0, v is a corresponding first-order KKT for θ. We see that because θl ⊙ β̂l = 0, βˆli = 0 when θli ̸= 0. This implies that we just need to show βˆli ≥ 0 when θli = 0. The second-order KKT criticality on smooth manifold implies that ∀ ∆θl s.t. 1T (ψ ′ (θl ) ⊙ ∆θl ) = 0, l = 1, . . . , L.   Consider tangent vectors of the type ∆θ = 0 . . . δl . . . 0 . Here δli is set to 1 only if θli = 0. All other components of δl including Euclidean tangent vectors are set to 0. Using this special tangent vector in the equation for the Product Hessian, we get ∆θ∇2 L2 ∆θ ≥ 0,

   0  ..  .. . .  i     2 ′′ s i    ∆θ∇ L2 ∆θ = 0 . . . δl . . . 0  diag(ψ (θl ) ⊙ (∇vl f + λl 1))   δl   . .  ..  . . 

0 ′

T

+(∆θ ⊙ φ (θ))

∇2v f (∆θ ⊙ φ′ (θ)). (25)

Noting that diag(ψ ′ (θl ))δl = 0 in view of the choice of δl , and recalling (4), we get ∆θT ∇2 L2 ∆θ = δlT diag (ψ ′′ (θl ) ⊙ ∇vl f + λsl ψ ′′ (θl )) δl = δlT diag (ψ ′′ (θl ) ⊙ (∇vl f + λsl 1)) δl   = δlT diag ψ ′′ (θl ) ⊙ β̂l δl = ψ ′′ (θl )i β̂li .

(26)

Title Suppressed Due to Excessive Length

11

Since the LHS is nonnegative by virtue of being a second-order KKT point and ψ ′′ (θl ) ≥ 0 because of our choice of convex parameterization, β̂li ≥ 0, l = 1, . . . , L. Also, ∇θq L2 = 0 =⇒ ∇vq L1 = 0. (27) Thus second-order criticality on smooth manifold implies first-order criticality on Fsimplex . ⊔ ⊓ We now strengthen Theorem 4.3 to get the converse of Theorem 4.2. Theorem 4.4 If θ⋆ is a second-order KKT point of (P×smooth ) then v ⋆ = φ(θ⋆ ) is a weak second-order KKT point of (Psimplex ). Proof Finally, we use our Hessian formula in (10) along with Li et al. [21, Equations 32-39] to prove that second-order critical point on smooth manifold implies weak second-order criticality on Fsimplex . ∆θli ψ ′ (θl ) = ∆vli

if

θli ̸= 0,

∆θli = 0

if

θli = 0.

(28)

It is easy to verify that this satisfies the smooth tangent space condition ⟨∆θl , ψ ′ (θl )⟩ = 0, l = 1, . . . , L. Using first-order criticality and the expression for Hessian on smooth manifold we get ∆θT ∇2v L2 ∆θ =

L X

∆θlT diag(β̂l ⊙ ψ ′′ (θl ))∆θl

l=1 T

+ ∆θ diag(φ′ (θ))∇2v f diag(φ′ (θ))∆θ =

L X ⟨∆θl , β̂l ⊙ ψ ′′ (θl ) ⊙ ∆θl ⟩ + ∆v T ∇2v f ∆v l=1

= 0 + ∆v T ∇2v f ∆v.

(29)

Here the first term becomes 0 since βˆli is 0 if θli is non zero, while ∆θli is 0 if θli is 0. Since the LHS is nonnegative, the RHS has to be nonnegative, proving our assertion. ⊓ ⊔ Remark 4.1 Theorem 4.4 becomes false if “weak” is replaced by “strong”. To see this, consider the simplex in R3 and the functions f (v) = −(v 1 − 1)2 and φ(θ) = θ ⊙ θ. The point θ = (1, 0, 0) is both a weak and strong second-order KKT point, but v = (1, 0, 0) is not a strong second-order KKT point. However, in view of Theorem 4.4, v is a weak second-order KKT point. In fact, regardless of the objective function, every vertex of the simplex is a weak second-order KKT point because, under condition (6), the constraints on ∆v imply ∆v = 0, and thus ∆v T ∇2 L1 ∆v = 0 ≥ 0 is satisfied.

12

Shashwat Kumar et al.

Proposition 4.1 Let ψ(θ)satisfy (1), and Im(ψ ′ ) = R, then its Legendre transform ψ ∗ (θ∗ ) = supθ∈R θ∗ θ − ψ(θ) defines an alternative reparameterization with the same second-order KKT–preserving properties. Proof Differentiating θ∗ θ − ψ(θ) with respect to θ and setting equal to zero we get θ∗ = ψ ′ (θ). Plugging this into the definition of Legendre transform we get ψ ∗ (ψ ′ (θ)) = ψ ′ (θ) θ − ψ(θ).

(30)

Setting θ = 0 and using ψ(0) = 0 yields ψ ∗ (ψ ′ (0)) = 0 and so ψ ∗ (0) = 0 because ψ ′ (0) = 0 according to remark 2.1. Differentiating ψ ∗ (θ∗ ) = θ∗ θ−ψ(θ) with respect to θ∗ gives (ψ ∗ )′ (θ∗ ) = θ and further differentiation gives (ψ ∗ )′′ (θ∗ ) =

dθ . dθ∗

(31)

Differentiating θ∗ = ψ ′ (θ) with respect to θ gives dθ∗ = ψ ′′ (θ), dθ

(32)

dθ 1 = ′′ . dθ∗ ψ (θ)

(33)

and hence by inversion

Therefore, (ψ ∗ )′′ (θ∗ ) =

1 ψ ′′ (θ)

> 0.

(34)

Therefore ψ ∗ satisfies the same conditions as in (1). As a consequence of Theorems 4.1, 4.2, 4.3, and 4.4, φ∗ (θ∗ ) defines a valid reparameterization that preserves second-order KKT properties. We showed above that ψ ∗ (0) = 0. Since ψ ∗ (θ∗ ) = supθ∈R {θ∗ θ − ψ(θ)}, define the Legendre transform in terms of θ as L(θ) := ψ ∗ (ψ ′ (θ)) = θψ ′ (θ) − ψ(θ). Since L′ (θ) = θψ ′′ (θ), and ψ ′′ (θ) > 0, L attains its unique minimum at θ = 0, where L(0) = 0. Also, limθ→±∞ L(θ) = ∞. Thus, Im(ψ ∗ ) = Im(L) = [0, ∞), so ψ ∗ is surjective onto [0, ∞), implying that the induced map to the simplex is surjective. ⊔ ⊓

5 Applications We now proceed to show two important applications of Theorems 4.1, 4.2, 4.3, and 4.4.

Title Suppressed Due to Excessive Length

13

5.1 Application 1: Simplex Structured Tensor Decomposition for Learning Probability Distributions Consider a set of N random variables {X1 , X2 , . . . , XN } where each Xn takes discrete values in [1, In ], n = 1, . . . , N . Let P represent the joint PMF tensor of these random variables where P(i1 , . . . , iN ) = Pr{X1 = i1 , . . . , XN = iN }. Similar to Kargas et al. [17], we assume that P admits a low-rank nonnegative Canonical Polyadic Decomposition (CPD) with F components, which means that there are N factor matrices A1 , . . . , AN and a loading vector λ such that: F X P= λϕ A1 (:, ϕ) ◦ A2 (:, ϕ) ◦ · · · ◦ AN (:, ϕ) ϕ=1

(35)

≜ [[λ, A1 , . . . , AN ]] 1T λ = 1, 1T An (:, ϕ) = 1,

n = 1, . . . , N,

ϕ = 1, . . . , F.

The problem is to fit model (35) to an empirical tensor Q estimated from data, according to some fit quality measure. This model can then be used to perform probabilistic topic modelling of sequences [25, 17] or estimate mutual information for doing feature selection [4]. Some previous works have focused on Least Squares Error [17] or Kullback-Leibler (KL) Divergence [31], and in our current work, we use the latter. We stack our factor matrices into a parameter vector v:
v = λ, A1 (:, 1), . . . , A1 (:, F ), . . . , AN (:, 1), . . . , AN (:, F ) , (36) v ∈ ∆F −1 × ∆I1 −1 × · · · × ∆I1 −1 × · · · × ∆IN −1 × · · · × ∆IN −1 . The KL-based objective function then reads f (v) =

IN X iN =1

···

I1 X i1 =1

Q(i1 , . . . , iN ) ln

Q(i1 , . . . , iN ) P̂(i1 , . . . , iN ; v)

(37)

where P̂(; v) indicates the tensor generated by unstacking the parameter v into the factor matrices and use the outer product in (35) to generate the PMF tensor P. Based on these operations, we introduce the RGD Algorithm, as shown in Algorithm 1 and RGD with Barzilai-Borwein step (RGD BB), outlined in Algorithm 2 for both applications of this paper. Algorithms 1 and 2 also require us to define product retraction and product transport. On the sphere Sn−1 , we define them by projection, respectively onto Sn−1 and onto its tangent space: retractionθ : Tθ Sn−1 → Sn−1 , θ + ∆θ retractionθ (∆θ) = , ∥θ + ∆θ∥ 2 product retraction = retraction × · · · × retraction

(38) (39)

14

Shashwat Kumar et al.

transportθ1 : Tθ Sn−1 → Tθ2 Sn−1 , transportθ1 (∆θ) = ∆θ − ⟨∆θ, θ2 ⟩ θ2 , product transport = transport × · · · × transport

(40) (41)

They are bonafide retraction and vector transport in view of [3, Proposition 5] and [2, (8.10)]. Algorithm 1 RGD with Backtracking Update Step Require: initial step size η > 0, initial reparameterized factor vector θ, empirical tensor Q, Armijo constant c > 0, and backtracking decay factor β ∈ (0, 1). 1: α ← η Q ▷ Use η as initial step size 2: GR ← θ ∇θ f (θ ⊙ θ) ▷ Riemannian Gradient via projection in product tangent space 3: m ← 1 4: ArmijoFlag ← False 5: while ArmijoFlag = False and m ≤ 25 do ▷ See, e.g., (39) 6: θ̂ = product retraction(θ, −α ∗ GR ) 7: if f (θ̂) − f (θ) ≤ −c ∗ α ∗ ∥GR ∥2 then 8: ArmijoFlag = True 9: else 10: α←β×α 11: end if 12: m←m+1 13: end while 14: return θ̂

5.2 Application 2: SRVF Registration for Functional Data Analysis In this section, we briefly discuss the SRVF registration framework, referring to [19] for further details. In contrast to [19], where SRVF registration is performed using dynamic programming and used for Shape PCA and clinical motion analysis, our goal here is to show that the registration problem has a product-simplex structure and can be optimized using the smooth reparameterization framework developed in this paper. As mentioned in [19], let {βl }L l=1 be a collection of curves representing some functional observations, each βl : [0, 1] → R. Let Γ denote the set of smooth, strictly increasing maps γ : [0, 1] → [0, 1] satisfying γ(0) = 0, γ(1) = 1, and γ̇ > 0. Elements of Γ are positive diffeomorphisms, also called phases, and are used to model temporal alignments. For a curve βl and a warping function γl ∈ Γ, the composition βl ◦ γl represents a time-warped version of βl . We first describe the alignment problem in the pairwise setting. Given two curves, β1 and β2 , the goal is to find a time-warping function γ2 such that the peaks and valleys of β2 ◦ γ2 are aligned with those of β1 . One way q to formulate R1 ⋆ f (t)2 dt this problem is γ2 ∈ arg minγ∈Γ ∥β1 − β2 ◦ γ∥2 where ∥f ∥2 = 0 2 denotes the classical L norm. However, this direct formulation has several

Title Suppressed Due to Excessive Length

15

Algorithm 2 RGD BB with Backtracking Update Step Require: initial step size η > 0, initial reparameterized factor vector θ, empirical tensor Q, Armijo constant c > 0, backtracking decay factor β ∈ (0, 1), Barzilai–Borwein step-size bounds αmin , αmax satisfying 0 < αmin < αmax < ∞, and memory parameter M ∈ N. 1: α ← η Q ▷ Use η as initial step size 2: GR ← θ ∇θ f (θ ⊙ θ) ▷ Riemannian Gradient via projection in product tangent space 3: m ← 1 4: ArmijoFlag ← False 5: while ArmijoFlag = False and m ≤ 25 do ▷ See, e.g., (39) 6: θ̂ = product retraction(θ, −α · GR ) 7: if f (θ̂) ≤ max f (θi ) − c ∗ α ∗ ∥GR ∥2 then i=max(1,m−M ),...,m

8: ArmijoFlag ← True 9: else 10: α←β·α 11: end if 12: m←m+1 13: end while Q 14: Gnew ← θ ∇θ f (θ̂ ⊙ θ̂) R 15: G̃R ← product transport(GR , θ, θ̂) 16: S ← −α · G̃R 17: Y ← Gnew − G̃R R 18: denom ← |⟨S, Y ⟩| 19: if denom < 10−30 then 20: αBB ← αmax 21: else ⟨S,S⟩ 22: αBB ← denom 23: end if 24: α ← max(min(αBB , αmax ), αmin ) 25: θ ← θ̂ 26: return θ̂, α

▷ New gradient ▷ See, e.g., (41)

▷ Return updated parameters and new step size

mathematical and practical limitations. In particular, this optimization is degenerate because ∥ · ∥2 is not invariant to reparameterization of time and ∥β1 − β2 ∥2 ̸= ∥β1 ◦ γ − β2 ◦ γ∥2 . The SRVF framework addresses this limitation by using an elastic Riemannian metric that is invariant to simultaneous time warpings of the input curves. The use of this metric is simplified using a mathematical representation called the square-root velocity function (SRVF). For a curve βl (t), the q

SRVF is defined as ql (t) = sign(β̇l (t))

|β̇l (t)| = √β̇l (t) . In case, β̇l (t) = 0, |β̇l (t)|

we set ql (t) = 0. Under a time warping of βl by γl , i.e., βl 7→ βl ◦ γl , the SRVF √ transforms as γl · ql = (ql ◦ γl ) γ̇l . This transformation is the basis of the elastic approach to curve alignment. Given two curves β1 and β2 , the pairwise SRVF registration problem is p 2 γ2⋆ = arg min q1 − (q2 ◦ γ) γ̇ . γ∈Γ

(42)

2

The corresponding SRVF distance is invariant to time warpings of the input curves, so it compares shapes rather than the rates at which the curves are traversed.

16

Shashwat Kumar et al.

For multiple curves, the registration problem can be written as ! L X p 2 µ̂L = min min µ − (ql ◦ γl ) γ̇l , µ

l=1

γl

(43)

where µ̂L is the mean SRVF shape and optimal γl aligns the lth curve to this mean. We now use this formulation to obtain an optimization problem over products of simplices. Let st = t/T , t = 0, . . . , T , be a uniform grid on [0, 1] with spacing ∆t = 1/T . Thus, the grid has T + 1 points and T intervals. For each warping function, define the discrete warping-rate increment variable vl (t) = γl (st+1 ) − γl (st ),

t = 0, . . . , T − 1.

(44)

Equivalently, when the warping function is represented through its derivative, vl (t) ≈ γ̇l (st )∆t. (45) Recall the boundary conditions γl , γl (0) = 0, γl (1) = 1, and the strictly increasing constraint γ̇l > 0. In practice, the strictly increasing constraint is relaxed with a non-decreasing γ̇l ≥ 0. This allows us to rewrite the boundary conditions in terms of the derivative: Z 1 γ̇l dt = γl (1) − γl (0) = 1. (46) t=0

For its discrete analog, this is equivalent to: T −1 X

vl (t) = 1.

(47)

t=0

Thus, each vl lies in the simplex ( T

∆ =

v ∈ RT≥0 :

T −1 X

) v(t) = 1 .

(48)

t=0

The corresponding discrete warping function is obtained by cumulative summation: γld (s0 ) = 0,

γld (st ) =

t−1 X

vl (k),

t = 1, . . . , T.

(49)

k=0

We represent the discretized mean SRVF by a Euclidean variable vL+1 ∈ RT . Hence, the feasible set for the discretized registration problem is Fsimplex = ∆T × · · · × ∆T ×RT , {z } |

(50)

L times

where the first L factors correspond to the warping-rate variables and the final factor corresponds to the discretized mean SRVF.

Title Suppressed Due to Excessive Length

17

Therefore, Theorems 4.1, 4.2, 4.3, and 4.4 apply to this discretized SRVF registration problem and we compare PGD, RGD, PGD BB, and RGD BB for optimizing the resulting registration objective. For each algorithm, we update vl in each iteration, and subsequently determine the corresponding warping function γld . By composing an unaligned function with corresponding γld using linear interpolation, we obtain the aligned function. The aligned functions are then used to update the mean vL+1 in each iteration. Through this iterative process, we ultimately achieve the optimal vl and the corresponding aligned functions. Autograd in PyTorch is used to obtained gradient of the objective function.

6 Numerical Results We benchmarked the following algorithms on a MacBook Air (Apple M3, 25GB RAM). All algorithms were implemented in Jax and PyTorch library in Python. 1. RGD BB with backtracking update step [16], outlined in Algorithm 2. 2. RGD with Armijo Condition based backtracking (Algorithm 1) as outlined in Li et al. [21], without the Wolfe condition. 3. PGD with projection based on Duchi et al. [10] outlined in Algorithm 3. 4. PGD with Barzilai-Borwein step size (PGD BB) [7] with projection based on Duchi et al. [10].

Algorithm 3 PGD with Backtracking Update Step Require: initial step size η > 0, initial reparameterized factor vector θ, empirical tensor Q, Armijo constant c > 0, and backtracking decay factor β ∈ (0, 1). 1: v = θ ⊙ θ 2: GE ← ∇v f (v) ▷ Euclidean Gradient 3: v̂ = product simplex projection(v − η ∗ GE ) 4: α ← 1 5: m ← 1 6: ArmijoFlag ← False 7: while ArmijoFlag = False and m ≤ 25 do 8: vnew = v + α ∗ (v̂ − v) 9: if f (vnew ) − f (v) <= c ∗ α ∗ ⟨GE , v̂ − v⟩ then 10: ArmijoFlag = True 11: else 12: α←β×α 13: end if 14: m←m+1 15: end while 16: return vnew

18

Shashwat Kumar et al.

Fig. 2 For each algorithm, we choose best step size η for different rank and tensor size. RGD BB outperforms other algorithm, often by several orders of magnitude. Higher-dimensional tensors are harder to factorize than lower-dimensional ones, although rank of tensor does not seem to affect the final quality of the solution as much.

6.1 Simplex Constrained Canonical Polyadic Tensor Factorization In this section, we present results from numerical simulations to demonstrate the performance of our algorithms. We generate a random joint probability distribution with three random variables, P (X1 , X2 , X3 ). The data is generated from a Naive Bayes CPD Model: P (X1 , X2 , X3 ) =

F X

P (X1 |t)P (X2 |t)P (X3 |t)P (t).

ϕ=1

The number of values of each variable Xn is chosen from the set {4, 8, 16}. The rank F of the CPD is selected from {4, 8, 16}. For all algorithms, we test initial step sizes (η) of the form 2i , where i ranges between [−4, 4]. All algorithms are run for a combination of these parameters until a maximum time of 20 seconds is reached. The decay factor β is set to 0.5 for both RGD BB and PGD BB, and 0.75 for RGD and PGD. A small decay constant c = 10−5 is chosen for the BB algorithms, while c = 10−4 is used for PGD and RGD. These parameters were hand-tuned to achieve optimal performance for each algorithm. The value m ≤ 25 in line 5 of Algorithm 1 comes from [21]. In Fig. 2, we report the results of the best performing initial step size η for each algorithm across a combination of tensor sizes and ranks. As seen in most rank and tensor size combinations, RGD BB performs the best, often

Title Suppressed Due to Excessive Length

19

outperforming the other algorithms by several orders of magnitude in the final solution cost. This improvement appears to be primarily due to RGD BB requiring fewer iterations to reach low-cost solutions. Interestingly, while higher-dimensional tensors were harder to factorize, as evidenced by the lower KL-Divergence, we did not observe a similar effect with the rank. In our Appendix A, we report results on different initial step sizes η for each algorithm. Out of all the algorithms, RGD BB often performs the best across a range of step sizes, being less sensitive to its misspecification. This makes it particularly attractive as a plug-and-play algorithm for large-scale problems, where extensive grid search over parameters is infeasible.

6.2 Functional Data Registration We start by randomly generating 10 sinusoidal functions βl with varying amplitudes (Fig. 3(a), first column). These could, for example, indicate a population of L individuals performing a periodic motion (e.g., doing a bicep curl) at different speeds. These amplitude functions are then composed with some random diffeomorphic warpings of time γ. The result is unaligned functions in the third column of Fig. 3(a), which can synthesize the effect of population level variability across individuals performing the same motion at different rates. Our goal is to recover the original sinusoidal shape from these temporally distorted motions. We parameterize the problem on the simplex in terms of v and also reparameterize it smoothly in terms of θ under ψ = θ ⊙ θ. The warping functions vl are randomly intialized on the simplex ensuring vl ≥ 0 and 1T vl = 1. The Euclidean mean SRVF parameter is stored in vL+1 . For all algorithms, we alternate between the following two steps: (a) computing the optimal warping function γld via optimization of vl , and (b) composing the unaligned functions

Fig. 3 Comparison of algorithms in SRVF registration: a) Starting from unaligned functions in third column, RGD and RGD BB results in much better alignment and reconstruction of original sinousoidal function βˆl shapes whereas PGD and PGD BB results in a distorted waveform. b) SRVF registration error curve of RGD and RGD BB compared to PGD and PGD BB.

20

Shashwat Kumar et al.

with the warping functions to obtain the aligned functions and updating the SRVF mean vL+1 . In this manner, all algorithms are run for 60 seconds per trial, and the reported results are averaged over 50 trials. To determine the optimal learning rate, we try step sizes in 2i , i ranging between [−4, 4]. Ultimately, we selected 0.5 as our choice of step size as it yielded the best results. For all the algorithms, the decay factor β is set to 0.5 and decay constant is chosen as c = 10−4 . Similar to the tensor decomposition problem, these parameters were hand-tuned to achieve optimal performance for each algorithm. Fig. 3(a) illustrates the qualitative alignment results. The RGD and RGD BB warping functions are smoother and align more closely with the original compared to those produced by PGD and PGD BB (Fig. 3(a), fourth column). As a result, the RGD and RGD BB aligned functions remain closer to the original sine waves (Fig. 3(a), fifth column). In contrast, the PGD and PGD BB warping functions display irregularities and roughness, leading to aligned functions that deviate from the original shapes. Specifically, there are abrupt changes at the edges in PGD, and as a result, the sine waves transform into square waves. This is likely due to the projection operation in PGD resulting in γ̇ ending up with solutions with zeros, where the reconstructed γ look much more step like. However, in the reparameterized problem, the natural smoothness of sphere results in smoother warping functions, which are critical in applications requiring accurate diffeomorphic shape registration. Fig. 3(b) presents the mean log registration error curves for all the approaches, where the error corresponds to the value of the SRVF registration objective function. It demonstrates that RGD and RGD BB achieve faster convergence and more effective optimization, resulting in lower error compared to PGD and PGD BB. These findings emphasize the advantages of Riemannian techniques in achieving efficient and reliable performance in aligning functional data.

7 Conclusion In this work, we presented a framework for smooth reparameterization of L vectors on simplicial product spaces. We showed that reparameterization preserves weak second-order KKT points and leveraged an efficient Riemannian Barzilai-Borwein algorithm to address the reparameterized problem. The Riemannian algorithm achieves excellent quantitative performance on both problems mentioned in section 6 and is robust to the choice of initial step size.

Funding This work was supported by the Fonds de la Recherche Scientifique – FNRS under Grant no T.0001.23.

Title Suppressed Due to Excessive Length

21

Competing Interests The authors declare that they have no competing interests.

Code Availability The implementation used to generate all numerical experiments in this paper is available at https://github.com/Arafat245/smooth_reparameterizations.

References 1. Abdolali, M., Gillis, N.: Simplex-structured matrix factorization: Sparsity-based identifiability and provably correct algorithms. SIAM Journal on Mathematics of Data Science 3(2), 593–623 (2021) 2. Absil, P.A., Mahony, R., Sepulchre, R.: Optimization Algorithms on Matrix Manifolds. Princeton University Press, Princeton, NJ (2008). URL http://press.princeton.edu/ titles/8586.html 3. Absil, P.A., Malick, J.: Projection-like retractions on matrix manifolds. SIAM Journal on Optimization 22(1), 135–158 (2012). DOI 10.1137/100802529. URL http://link. aip.org/link/?SJE/22/135/1 4. Amiridi, M., Kargas, N., Sidiropoulos, N.D.: Information-theoretic feature selection via tensor decomposition and submodularity. IEEE Transactions on Signal Processing 69, 6195–6205 (2021) 5. Aragón, F.J., et al.: Nonlinear Optimization. Society for Industrial and Applied Mathematics (2019) 6. Bertsekas, D.P.: Nonlinear programming. Journal of the Operational Research Society 48(3), 334–334 (1997) 7. Birgin, E.G., Martı́nez, J.M., Raydan, M.: Nonmonotone spectral projected gradient methods on convex sets. SIAM Journal on Optimization 10(4), 1196–1211 (2000) 8. Bonnans, J.F., Gilbert, J.C., Lemaréchal, C., Sagastizábal, C.A.: Numerical optimization: theoretical and practical aspects. Springer Science & Business Media (2006) 9. Boumal, N.: An introduction to optimization on smooth manifolds. Cambridge University Press (2023). DOI 10.1017/9781009166164. URL https://www.nicolasboumal. net/book 10. Duchi, J., Shalev-Shwartz, S., Singer, Y., Chandra, T.: Efficient projections onto the l1ball for learning in high dimensions. In: Proceedings of the 25th international conference on Machine learning, pp. 272–279 (2008) 11. Esposito, F., Ang, A.: Riemannian multiplicative update for sparse simplex constraint using oblique rotation manifold. arXiv preprint arXiv:2503.24075 (2025) 12. Faybusovich, L.: Hamiltonian structure of dynamical systems which solve linear programming problems. Physica D: Nonlinear Phenomena 53(2-4), 217–232 (1991) 13. Gill, P.E., Murray, W., Wright, M.H.: Practical optimization. SIAM (2019) 14. Guo, L., Lin, G.H., Ye, J.J.: Second-order optimality conditions for mathematical programs with equilibrium constraints. Journal of Optimization Theory and Applications 158, 33–64 (2013). DOI 10.1007/s10957-012-0228-x 15. Guo, Z., Min, A., Yang, B., Chen, J., Li, H., Gao, J.: A sparse oblique-manifold nonnegative matrix factorization for hyperspectral unmixing. IEEE Transactions on Geoscience and Remote Sensing 60, 1–13 (2021) 16. Iannazzo, B., Porcelli, M.: The Riemannian Barzilai–Borwein method with nonmonotone line search and the matrix geometric mean computation. IMA Journal of Numerical Analysis 38(1), 495–517 (2018)

22

Shashwat Kumar et al.

17. Kargas, N., Sidiropoulos, N.D., Fu, X.: Tensors, learning, and “Kolmogorov extension” for finite-alphabet random vectors. IEEE Transactions on Signal Processing 66(18), 4854–4868 (2018) 18. Kizel, F., Shoshany, M., Netanyahu, N.S., Even-Tzur, G., Benediktsson, J.A.: A stepwise analytical projected gradient descent search for hyperspectral unmixing and its code vectorization. IEEE Transactions on Geoscience and Remote Sensing 55(9), 4925–4943 (2017) 19. Kumar, S., Rahman, A., Gutierrez, R., Livermon, S., McCrady, A.N., Blemker, S., Scharf, R., Srivastava, A., Barnes, L.E.: A shape-based functional index for objective assessment of pediatric motor function. Plos one 20(10), e0332,383 (2025) 20. Levin, E., Kileel, J., Boumal, N.: The effect of smooth parametrizations on nonconvex optimization landscapes. arXiv preprint arXiv:2207.03512 (2022) 21. Li, Q., McKenzie, D., Yin, W.: From the simplex to the sphere: Faster constrained optimization using the Hadamard parametrization. arXiv preprint arXiv:2112.05273 (2021) 22. Ma, W.K., Bioucas-Dias, J.M., Chan, T.H., Gillis, N., Gader, P., Plaza, A.J., Ambikapathi, A., Chi, C.Y.: A signal processing perspective on hyperspectral unmixing: Insights from remote sensing. IEEE Signal Processing Magazine 31(1), 67–81 (2013) 23. Niculescu, C., Persson, L.E.: Convex functions and their applications, vol. 23. Springer (2006) 24. Nocedal, J., Wright, S.J.: Numerical Optimization, 2nd edn. Springer (2006) 25. Sidiropoulos, N.D., De Lathauwer, L., Fu, X., Huang, K., Papalexakis, E.E., Faloutsos, C.: Tensor decomposition for signal processing and machine learning. IEEE Transactions on signal processing 65(13), 3551–3582 (2017) 26. Sidiropoulos, N.D., Karakasis, P.A., Konar, A.: Minimizing low-rank models of highorder tensors: Hardness, span, tight relaxation, and applications. IEEE Transactions on Signal Processing (2023) 27. Srivastava, A., Klassen, E.P.: Functional and shape data analysis, vol. 1. Springer (2016) 28. Tucker, J.D., Wu, W., Srivastava, A.: Generative models for functional data using phase and amplitude separation. Computational Statistics & Data Analysis 61, 50–66 (2013) 29. Vaskevicius, T., Kanade, V., Rebeschini, P.: Implicit regularization for optimal sparse recovery. Advances in Neural Information Processing Systems 32 (2019) 30. Wu, R., Ma, W.K., Fu, X.: A stochastic maximum-likelihood framework for simplex structured matrix factorization. In: 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 2557–2561. IEEE (2017) 31. Yeredor, A., Haardt, M.: Estimation of a low-rank probability-tensor from sample subtensors via joint factorization minimizing the Kullback-Leibler divergence. In: 2019 27th European Signal Processing Conference (EUSIPCO), pp. 1–5. IEEE (2019) 32. Zhao, P., Yang, Y., He, Q.C.: Implicit regularization via Hadamard product overparametrization in high-dimensional linear regression. arXiv preprint arXiv:1903.09367 2(4), 8 (2019)

Title Suppressed Due to Excessive Length

23

A Appendix (Results with Different Learning Rates)

Fig. 4 Results for tensor dimensions 4, 8, 16, 32 of rank 4 with different learning rates η.

24

Shashwat Kumar et al.

Fig. 5 Results for tensor dimensions 4, 8, 16, 32 of rank 8 with different learning rates η.

Fig. 6 Results for tensor dimensions 4, 8, 16, 32 of rank 16 with different learning rates η.

Related documents

• Borel Isomorphism of SPR Markov Shifts#395228
• Keep the Angle: A Universal Geometry-Preserving Basis in Spectral Embeddings#390556
• A multi-dimensional quantum estimation and model learning framework based on variational Bayesian inference#149325
• Tensor Decomposition for Energy-Momentum Correlation Functions#1109075
• Determination of all diffeomorphism invariant tensor fields on the space of smooth positive densities on a compact manifold with corners#1128864