Generalized cyclic symmetric decompositions for the matrix multiplication tensor
Charlotte Vermeylen Marc Van Barel
































label1]Charlotte Vermeylen
[label1]organization=Department of Electrical Engineering (ESAT) KU Leuven,
             addressline=Kasteelpark Arenberg 10,
             city=Heverlee,
             postcode=3001,
             state=Vlaams-Brabant,
             country=Belgium
             
label2]Marc Van Barel
[label2]organization=Department of Computer Science KU Leuven,
             addressline=Celestijnenlaan 200C A,
             city=Heverlee,
             postcode=3001,
             state=Vlaams-Brabant,
             country=Belgium











In this paper we propose a new generalized cyclic symmetric structure in the factor matrices of polyadic decompositions of matrix multiplication tensors for non-square matrix multiplication to reduce the number of variables in the optimization problem and in this way improve the convergence.















matrix multiplication tensors optimization canonical polyadic decomposition














§ INTRODUCTION
 
In this paper we consider the fast matrix multiplication (FMM) problem which rewrites matrix multiplication, which is a set of bilinear equations, i.e.,

    C(i,j) = ∑_k=1^p A(i,k) B(k,n),

where C:=AB, as a tensor equation:

    C    = ∑_i,j,k=1^m,p,n⟨ E_ik^m× p, A  ⟩_F⟨ E_kj^p× n, B ⟩_F E_ij^m× n
       = ( ∑_i,j,k=1^m,p,n E_ik^m× p⊗ E_kj^p× n⊗ E_ij^m× n)·_F( A, B ),

where E_ij^m× p, for i=1,…,m and j=1,…, p, is a basis of matrices in ℝ^m× p such that E_ij^m× p(i_1,i_2) = 1 if (i_1,i_2)=(i,j) and zero else, `⊗' denotes the tensor product and ⟨·⟩_F and `·_F' denote the Frobenius inner product. We add a transpose to (<ref>) such that the tensor that is defined has additional interesting properties, such as cyclic symmetry (CS), which will be discussed in more detail in <Ref>: 

    C^⊤   = ( ∑_i,j,k=1^m,p,n E_ik^m× p⊗ E_kj^p× n⊗(E_ij^m× n)^⊤)·_F( A, B ) 
       = ( ∑_i,j,k=1^m,p,n E_ik^m× p⊗ E_kj^p× n⊗ E_ji^n× m)·_F( A, B ).

Remark that the Frobenius inner product is a linear function in the elements of a matrix. The rank of the bilinear equation is the minimal r  such that the equation can be written as

    C^⊤   =  ∑_i=1^r⟨ U_i, A  ⟩_F⟨ V_i, B ⟩_F W_i,

for some matrices U_i ∈ℝ^m × p, V_i ∈ℝ^p × n, and W_i ∈ℝ^n × m. It can be shown that minimizing r minimizes the computational complexity of matrix multiplication <cit.>, <cit.>. More specifically, the number of arithmetic operations needed to compute two matrices in ℝ^n × n is 𝒪( n^ω), where ω:= log_n r. Since r < n^3, it holds that ω < 3. Additionally, it is well known that the exponent is bounded from below by two, since you need at least n^2 operations to compute all n^2 elements of the matrix multiplication.

Minimizing the rank of the bilinear equation corresponds to finding a (canonical) polyadic decomposition ((C)PD) of the so called matrix multiplication tensor (MMT) <cit.> defined implicitly in (<ref>):

    T_mpn   :=  ∑_i,j,k=1^m,p,n e_ik^mp⊗ e_kj^pn⊗ e_ji^mn,

where e_ik^mp:= vec(E_ik^m× p) and T_mpn is the MMT of dimension mp × pn × mn. Consequently, (<ref>) can also be formulated as

    vec (C^⊤) = T_mpn·_1 vec(A) ·_2 vec(B),

where vec( ·) denotes the column-wise vectorization operator, and ·_1 and ·_2 denote the multiplication along the first and second dimension of T_mpn respectively.

Because of this definition, T_mpn is a sparse tensor consisting of mpn ones: 

    T_mpn(i_1 + (i_2-1)m,j_1 + (j_2-1)p,k_1 + (k_2-1)n) 
       =  1     if   i_1 = k_2, i_2 = j_1,j_2 = k_1 
    
    0    else,

for all i_1,k_2 = 1, …, m, i_2,j_1 = 1, …, p, and j_2, k_1 = 1, …, n.
Remark that a MMT only depends of the size of the matrices that are multiplied and not on the values. 

A PD of T_mpn of rank r is a decomposition of T_mpn into r rank-1 tensors:

    T_mpn = ∑_i=1^r u_i ⊗ v_i ⊗ w_i,

where `⊗' denotes the tensor product and u_i,v_i, and w_i are vectors in ℝ^n_1, ℝ^n_2, and ℝ^n_3 respectively. These vectors are collected in three factor matrices U, V, and W:
U := [ u_1, ⋯, u_r ], V := [ v_1, ⋯, v_r ], and W := [ w_1, ⋯, w_r]. 
A PD of rank r is written in short as PD_r and a PD of rank r of a tensor T as PD_r(T). The minimal r for which a PD of a certain tensor exists is called the rank r^* of this tensor or r^*(T) and a PD of this rank a CPD or PD_r^*(T). 

With a PD_r (<ref>) of T_mpn, we obtain the following base algorithm to multiply any two matrices A ∈ℝ^m× p and B ∈ℝ^p × n:

    vec( C^⊤)    = ∑_i=1^r⟨ u_i, vec(A) ⟩⟨ v_i,  vec(B) ⟩ w_i,

which is the same as (<ref>) when U_i = reshape(u_i,[m,p]), V_i = reshape(v_i,[p,n]), and W_i = reshape(w_i,[n,m]), for all i=1, …, r. The number of multiplications between (linear combinations of) elements of A and B in (<ref>)  is reduced to r < mpn, where the upper bound holds for standard matrix multiplication. Such multiplications are called ‘active' because, when these algorithms are applied recursively, they determine the asymptotic complexity of matrix multiplication <cit.>.

One of the difficulties is that the rank of a CPD of T_mpn, denoted by r^*(T_mpn), is not known for most combinations of (m,p,n). Counterexamples are T_222, for which the set of PDs is completely understood and r^* is known to be seven <cit.>, T_223, for which the rank is known to be 11 <cit.>, and T_224, for which the rank is known to be 14 <cit.>. We call PD_7(T_222) originally discovered by Strassen PD_Strassen <cit.>. For other combinations of m, p, and n, lower bounds on the rank exist but none of these lower bounds are attained in practice. For example, the lower bound on r^*(T_333) is proven to be 19 <cit.>, and the upper bound is 23 <cit.>. We call the lowest rank for which an exact PD of T_mpn is known in the literature r̃(T_mpn). For an overview of r̃(T_mpn) for various m, p, and n, we refer to <cit.> and <cit.>. Remark that the rank of T_mpn is equal to the rank of T_pnm or any other permutation of m, p, and n. This is true because we can obtain a PD_r of T_pnm and T_nmp by cyclically permuting the factors in the PD:

    T_pnm = ∑_i=1^r v_i ⊗ w_i ⊗ u_i,        T_nmp = ∑_i=1^r w_i ⊗ u_i ⊗ v_i,

and additionally, it can be shown that
T_mnp = ∑_i=1^rvec(W_i^⊤) ⊗vec(V_i^⊤) ⊗vec(U_i^⊤), and similarly for T_npm and T_pmn.

To find a PD_r of T_mpn, we use the following nonlinear least squares (NLS) cost function, which is also frequently used in the literature for the FMM problem <cit.>:

    min_x  1/2‖  F(x)- vec(T_mpn) ‖^2_=: f(x),

where ‖·‖ denotes the ℓ_2-norm, and F(x) is a vector function defined as

    F(x )    :=  vec( ∑_r=1^r u_r ⊗ v_r ⊗ w_r ),     x    := [ vec(U); vec(V); vec(W) ].

Remark that f(x^*)=0, if and only if x^* is a PD_r of T_mpn. However, most standard optimization algorithms fail to converge to a global minimum of (<ref>). One of the reasons is that PDs of MMT have additional continuous invariances compared to generic PDs, which we will review in <Ref> <cit.>.

In practice, the minimal length or rank r^* is approximated experimentally by the lowest length r̃ for which a solution x^* with f(x^*)=0 can be obtained using numerical optimization. However, no globally convergent algorithm for (<ref>) is known to us, and thus we can never be sure that no solution with a lower rank exists. 

Since the total number of elements of the factor matrices grows rapidly with the size of the matrices, solving (<ref>) is only feasible for relatively small values of m, p, and n, e.g., m,p,n ≤ 5.
In this paper we give a new structure that can be enforced in PDs of MMTs to reduce the number of variables in the optimization problem and hereby improve the convergence.

To solve (<ref>), in the literature the alternating least squares (ALS) method is frequently used <cit.> but the convergence is usually slow and a lot of starting points are needed to obtain a solution in a reasonable amount of time. In <cit.> and <cit.>, a constrained optimization problem and corresponding method are proposed to improve the convergence. More specifically, a Levenberg-Marquardt (LM) method with Lagrange parameters and a quadratic penalty (QP) term in combination with the ALS method are used respectively. However, both methods were still not able to find PDs of rank 49 of MMTs for 4 × 4 matrix multiplication whereas these PDs are known to exist. That is why we proposed in <cit.> an augmented Lagrangian (AL) method and a new constrained optimization problem to find PDs of MMTs:

    min_x f(x),   s.t.   h(x)=0,    l ≤ x ≤ u.

Different equality constraints h(x) were proposed.
The AL method can be considered as a combination of the LM and QP method since the AL objective function is:

    min_x,y_1,y_2,zf(x) + ⟨ y_1,h(x) ⟩ + ⟨ y_2, x - z ⟩ + β/2( ‖ h(x) ‖^2 + ‖ x- z‖^2 )_=:ℒ_A ( x,y_1, y_2,z;β),

subject to (s.t.) l ≤ z ≤ u, where y_1 and y_2 are vectors of Lagrange multipliers, β is the regularization parameter, and z is a vector of slack variables that enables us to rewrite the inequality constraint as an equality constraint. Because f(x) is an NLS function, the AL objective can also be rewritten as an NLS function. The variables in x and the auxiliary parameters are updated successively according to <cit.>. We use the LM method for the minimization to x. The pseudo-code of the algorithm is shown in <cit.> and the  implementation is publicly available [<https://github.com/CharlotteVermeylen/ALM_FMM_stability>].
The advantage compared to the QP method is that the constraints are satisfied more accurately during optimization. Using this method, we were able to obtain new PD_49(T_444)s and also different new CS PDs. We will use this method in the numerical experiments in <Ref>.

The rest of the paper is organized as follows. In <Ref>, we give some preliminaries concerning the FMM problem. In <Ref>, the new generalized CS structure is proposed, and in <Ref>, numerical experiments are given to demonstrate the use and advantages of including this structure in the optimization problem (<ref>).





§ PRELIMINARIES
 

In this section, we first discuss what we mean with practical PDs and FMM algorithms. Afterwards, we give some well known properties of MMTs.



 §.§ Practical algorithms

As discussed in the introduction, a PD_r(T_mpn) can be applied recursively to, e.g., multiply matrices of size m^k × p^k and p^k × n^k. Such a recursive algorithm requires 𝒪( c r^k) operations, where the constant c depends on the number of scalar multiplications and nonzeros in the factor matrices. When k is sufficiently large and c sufficiently small, fewer operations are required compared to standard matrix multiplication. To decrease the constant c, we search for sparse PDs, with elements in a discrete set, e.g., {-1,0,1 }. This set can be extended with powers of 2 since this is not a costly operation when implemented in hardware. Such PDs are called practical PDs as they result in practical FMM algorithms. 

Consequently, if we obtain a numerical PD of T_mpn using numerical optimization, it still has to be transformed into a practical one. The invariance transformations or inv-transformations discovered in <cit.>, which are further discussed in <Ref>, can be used to obtain such PDs <cit.>. However, not all PDs of T_mpn can be transformed into a practical PD in this way <cit.>. In <cit.>, this process is called discretization.

A disadvantage of (<ref>) is that the solutions have floating point elements and thus extra steps have to be taken or constraints have to be added to the optimization problem to obtain practical PDs. That is why we proposed a new constrained problem formulation in <cit.> with equality constraint:

    h_discr (x):= x· (x-1) · (x+1),

where `·' denotes element-wise multiplication.



 §.§ Properties of the matrix multiplication tensor


In this section, first the invariance transformations that are present for all PDs of T_mpn are discussed. Additionally, a short overview of what is known in the literature about the use of these invariance transformations to obtain practical PDs is given. Afterwards, cyclic symmetry and recursive PDs are discussed in more detail. 



  §.§.§ Invariance transformations
 

The representation of a PD via factor matrices is not unique. For example when permuting the rank-1 tensors, and consequently the columns of the factor matrices according to

    ( U,V,W ) →( U',V',W' ),

where u'_i= u_σ(i), for all i=1,…, r, and σ = ( σ(1)    σ(2)    σ(3)    …   σ(r)
) is an element of the permutation group of r elements, then the same tensor is obtained:

    ∑_i=1^r u_i ⊗ v_i⊗ w_i = ∑_i=1^r u_i' ⊗ v_i'⊗ w_i'.

Additionally, because of the multi-linearity of the tensor product, the columns of each rank-1 tensor can be scaled as: α_i u_i, β_i v_i, 1/α_i β_i w_i, for all i=1, …, r, and α_i,β_i in ℝ_0, without changing the rank-1 tensors because

    α_i u_i ⊗β v_i  ⊗1/α_i β_i w_i = u_i ⊗  v_i  ⊗ w_i.

These two invariances hold for all PDs. Remark that only the scaling invariance is a continuous invariance transformation (of dimension 2r).

For PDs of T_mpn, additional invariances hold <cit.>. More specifically, it is well known that these PDs are invariant under the following PQR-transformation:

    u_i' ←vec( P U_i Q^-1),    v_i' ←vec( Q V_i R^-1),     w_i' ←vec( R W_i P^-1),

for i=1, …, r, and where P ∈ GL(m), Q ∈ GL(p), and R ∈ GL(n), where GL(i) represents the general linear group of invertible matrices in ℝ^i × i. This invariance can easily be proven using the fact that C := AB = (P^⊤)^-1( P^⊤ A (Q^⊤)^-1) ( Q^⊤ B (R^⊤)^-1) R^⊤. When we substitute this in (<ref>), we obtain:

    C = ∑_i=1^r    trace( U_i (P^⊤ A (Q^⊤)^-1)^⊤) 
       trace( V_i ( Q^⊤ B (R^⊤)^-1)^⊤) (P^⊤)^-1 W_i^⊤ R^⊤
    
    = ∑_i=1^r    trace( P U_i Q^-1 A^⊤) trace( Q V_i R^-1 B^⊤) ( R W_i P^-1)^⊤,

where we made use of the fact that the trace is invariant under cyclic permutation of the factors. And thus indeed another base algorithm or PD_r of T_mpn can be obtained using (<ref>). Additionally, when m=n=p, PDs of T_mmm are invariant under the following transpose-transformation:

    u_i' ←vec( V_i^⊤) ,    v_i' ←vec( U_i^⊤) ,    w_i' ←vec( W_i^⊤),

for i=1, …, r, which can be proven using the fact that C = ( B^⊤ A^⊤)^⊤. Lastly, still for m=n=p, the PDs are invariant under cyclic permutation of the factor matrices:

    ∑_i=1^r u_i ⊗ v_i⊗ w_i = ∑_i=1^r v_i ⊗ w_i⊗ u_i = ∑_i=1^r w_i ⊗ u_i⊗ v_i,

which holds because T_mmm is a CS tensor.

Remark that concerning the specific invariances of PDs of T_mpn, only the PQR-invariance is a continuous invariance and it overlaps with the scaling invariance when P=aI, Q=bI, and R=cI, where I is the identity matrix of the appropriate size and a,b, and c are scaling factors in ℝ_0, because in this case the PQR-transformation scales the columns of the factor matrix U with a/b, all columns of V with b/c, and all columns of W with c/a, and indeed the product of these factors equals one as for the scaling invariance. Thus, the combination of the scaling and PQR-transformation has at most dimension 2r + m^2 + p^2 + n^2 - 3 for all PDs of T_mpn.

We call the combination of the invariance transformations described above inv-transformations. Two PDs of T_mpn that can be obtained using inv-transformations are called inv-equivalent. Only for T_222, it is known that all PDs are inv-equivalent <cit.>.



  §.§.§ Cyclic symmetry
 

Because of its definition, T_mpn is a structured and more specifically cyclic symmetric (CS) tensor when m=p=n. This means that T_mmm(i,j,k) = T_mmm(j,k,i)=T_mmm(k,i,j), for all i,j,k=1,…, m^2. 

We can make use of this property for constructing a PD by requiring that the rank-1 tensors are symmetric or occur in CS pairs <cit.>:

    T_mmm = ∑_i=1^s a_i ⊗ a_i ⊗ a_i +∑_j=1^t ( b_j ⊗ d_j ⊗ c_j + c_j ⊗ b_j ⊗ d_j + d_j ⊗ c_j ⊗ b_j ),

where a_i,b_j,c_j, and d_j are vectors in ℝ^m^2, for all i =1,…, s and j=1,…, t, where s and t are parameters denoting respectively the number of symmetric and CS pairs of asymmetric rank-1 tensors. The length then equals r=s+3t. If we define the matrices
A := [ a_1⋯ a_s], 
B := [ b_1⋯ b_t], C := [ c_1⋯ c_t], and D := [ d_1⋯ d_t],
the factor matrices can be written in function of these matrices:

    U   = [ A B C D ],    V   = [ A D B C ],    W    = [ A C D B ].

Consequently, the number of unknowns is reduced by a factor 3, which makes this CS structure very useful for larger problem parameters, e.g., m,p,n≥ 4. Furthermore, the cost function can then be changed to

    min_A,B,C,D1/2(    ∑_i=1^m^2∑_j=i^m^2∑_k=i+1^m^2( T_mmm(i,j,k) - ∑_i'=1^s a_ii' a_ji' a_ki'
    
    -    ∑_j'=1^t ( b_ij' d_jj' c_kj' + c_ij' b_jj' d_kj' + d_ij' c_jj' b_kj') )^2 
    
    +    ∑_i=1^m^2( T_mmm(i,i,i) - ∑_i'=1^s a_ii'^3 - 3 ∑_j'=1^t ( b_ij' d_ij' c_ij') )^2 ),

where a_ii':=A(i,i'), b_ii':=B(i,i') and similarly for c_ii' and d_ii', which reduces the number of rows in the Jacobian matrix used in NLS optimization methods to solve (<ref>) from m^6 to 1/3 (m^6 - m^2) + m^2. Remark that this cost function only includes one element for each CS pair of points because the CS structure already ensures that these elements are equal.
A disadvantage of including this structure may be that by restricting the search space we might not be able to find the most sparse or stable PDs of T_mmm. Furthermore, it is not known if the rank of a PD of T_mmm with this structure equals the canonical rank.

Different practical CS PDs of rank 7 for T_222 and of rank 23 for T_333 are known in the literature <cit.>. These PDs were found using an ALS method and further investigated using algebraic geometry and group theory. Strassen's decomposition satisfies s=1 and t=2.
Although all PD_7( T_222)s are known to be inv-equivalent, PD_Strassen can be transformed into a practical PD with CS parameters (s,t)=(4,1). However, this PD contains more nonzeros than PD_Strassen and thus is not used in practice. 



 §.§ Decompositions obtained by recursion


As mentioned in the introduction, an FMM algorithm applies a base algorithm for small m, p, and n, to large matrices, e.g., of size m^k × p^k and p^k × n^k. In the same way, we can also obtain PDs of rank 49 of T_444 using a PD_7(T_222) one time recursively. The PD_49(T_444) that is obtained using PD_Strassen is called PD_Strassen^rec in the rest of the text. In the following proposition, we give the formulas to obtain the factor matrices of a recursive PD from the factor matrices of the original PD. This result is well known, e.g., from the supplementary material of <cit.>, and in <cit.>, and a proof can be found, e.g., in <cit.>.

 
If U, V, and W are the factor matrices of a PD_r(T_mpn), then the factor matrices U', V', and W', where

    u'_i'   := 2⃗( U_i_1⊗ U_i_2),     v'_i'   := 2⃗( V_i_1⊗ V_i_2), 
    
    w'_i'   := 2⃗( W_i_1⊗ W_i_2),     i'   :=i_2+(i_1-1)r,

for i_1,i_2 = 1, …, r, are the factor matrices of a PD_r^2(T_m^2p^2n^2).


The following corollary gives the CS parameters and factor matrices of a recursive PD as a function of the original PD. 
 
  If U, V, and W are CS factor matrices of a PD_r( T_mpn) with CS parameters s and t, then the factor matrices U', V', and W' of the recursive decomposition PD_r^2^rec( T_m^2p^2n^2) is also CS with parameters s'=s^2 and t'=t(s+r): 

    U'    := [ A'  B'  C'  D' ],     V'    := [ A'  D'  B' C' ],     W'    := [ A' C' D' B' ],

where 

    A'_i'   := A_i_1⊗ A_i_2,     i'    := i_2 + (i_1-1)s,     i_1,i_2    = 1, …, s,

and 

    B'    := [ B_1' B_2' ],     C'    := [ C_1' C_2' ],     D'    := [ D_1' D_2' ],

where 

    B'_1,j'   := A_i_1⊗ B_j,     C'_1,j'   := A_i_1⊗ C_j,     D'_1,j'   := A_i_1⊗ D_j, 
    
    B'_2,k'   := B_j⊗ U_k,     C'_2,k'   := C_j⊗ W_k,     D'_2,k'   := D_j⊗ V_k,

where j' := j + (i_1-1)t and k':=k + (j-1)r, for j = 1, …, t and k = 1 ,…, r. 


Using <Ref>, we know that U_i'' = U_i_1⊗ U_i_2, V_i'' = V_i_1⊗ V_i_2, and W_i'' = W_i_1⊗ W_i_2, where i':= i_2 + (i_1-1)r. Thus, the symmetric part of the recursive PD must satisfy: U_i_1 = V_i_1 = W_i_1 and U_i_2 = V_i_2  = W_i_2. Consequently, i_1 and i_2 must be smaller than or equal to s, such that U_i_1 = V_i_1 = W_i_1 = A_i_1 and U_i_2 = V_i_2  = W_i_2 = A_i_2. 
For the other values of i_1 and i_2, the columns still appear in CS pairs and they can be rearranged as in (<ref>) to satisfy the CS structure. In <cit.>, the different columns that are obtained by using <Ref> are written out to see this more clearly. 


Remark that, since PD_Strassen satisfies s=1 and t=2, the parameters of PD_Strassen^rec are s=1 and t=16. Another CS PD_49( T_444) can be obtained with (s,t) = (16,11) by using a PD_7( T_222) with (s,t) = (4,1) one time recursively.



 §.§ Jacobian matrix


In <cit.>, we showed that the Jacobian matrix can be used to investigate the inv-equivalence of PDs of MMTs and to investigate any additional invariances, such as the ones discovered in <cit.> using parametrizations of decompositions. That is why in <Ref>, we will give the size and ranks of the Jacobian matrix at solutions with the structure proposed in the next section. 



§ GENERALIZATION CYCLIC SYMMETRY
 

In this section, we give a generalization of the CS structure discussed in <Ref>, to the case when m, n, and p are not equal. We do this to reduce the number of variables and reduce the search space to hopefully speed up the convergence. Note however that similarly as for the CS structure discussed in <Ref>, it is not known if the rank of PDs with this structure equals the canonical rank. Because we can obtain a PD of T_pnm or another permutation of m, p, and n from a PD of T_mpn, we assume in this chapter that m ≤ p ≤ n. We know from (<ref>) that:

    T_mpn(i_1 + (i_2-1)m,j_1 + (j_2-1)p,k_1 + (k_2-1)n) 
       = ∑_i=1^r U_i(i_1,i_2) V_i(j_1,j_2) W_i(k_1,k_2) =  1     if   i_1 = k_2, i_2 = j_1,j_2 = k_1 
    
    0    else,

for all i_1,k_2=1,…, m, i_2,j_1=1,…, p, and j_2,k_1=1,…, n. 

Consequently, as long as i_1,i_2,j_1,j_2,k_1,k_2 ≤ m, it holds that

    T_mpn(i_1 + (i_2-1)m,j_1 + (j_2-1)p,k_1 + (k_2-1)n) 
    
    =    T_mpn(k_1 + (k_2-1)m,i_1 + (i_2-1)p,j_1 + (j_2-1)n) 
    
    =    T_mpn(j_1 + (j_2-1)m,k_1 + (k_2-1)p,k_1 + (k_2-1)n).

And thus,

    ∑_i=1^r U_i(i_1,i_2) V_i(j_1,j_2) W_i(k_1,k_2)    = ∑_i=1^r U_i(k_1,k_2) V_i(i_1,i_2) W_i(j_1,j_2)
       = ∑_i=1^r U_i(j_1,j_2) V_i(k_1,k_2) W_i(i_1,i_2),

which can be used to include a CS structure in a subpart of the factor matrices. An illustration is shown in <Ref>, where the following submatrices are defined in function of the position in the factor matrices:

    U_i(i_1,i_2) =: 
    A_i(i_1,i_2)    if  i_2 ≤ m, i ≤ s
    
    B_i-s(i_1,i_2)    if  i_2 ≤ m, s < i ≤ s+t
    
    C_i-s-t(i_1,i_2)    if  i_2 ≤ m, s+t < i ≤ s+2t
    
    D_i-s-2t(i_1,i_2)    if  i_2 ≤ m, s+2t < i ≤ r_CS
    Ũ_i(i_1,i_2-m)    if  i_2 > m,  i ≤ r_CS
    U̇_i-r_CS(i_1,i_2)    if  r_CS < i 
    ,

where r_CS:=s +3t, and for i=1,…, r, i_1=1,…, m, and i_2 = 1,…, p, 

    V_i(j_1,j_2) =: 
    A_i(j_1,j_2)    if  j_1,j_2 ≤ m, i ≤ s
    
    D_i-s(j_1,j_2)    if  j_1,j_2 ≤ m, s < i ≤ s+t
    
    B_i-s-t(j_1,j_2)    if  j_1,j_2 ≤ m, s+t < i ≤ s+2t
    
    C_i-s-2t(j_1,j_2)    if  j_1,j_2 ≤ m, s+2t < i ≤ r_CS
    V̂_i(j_1-m,j_2)    if  j_1 > m,j_2 ≤ m,  i ≤ r_CS
    Ṽ_i(j_1,j_2-m)    if  j_2 > m,  i ≤ r_CS
    V̇_i-r_CS(j_1,j_2)    if  r_CS < i 
    ,

for j_1=1,…, p, and j_2 = 1,…, n, and

    W_i(k_1,k_2) =: 
    A_i(k_1,k_2)    if  k_1 ≤ m, i ≤ s
    
    B_i-s(k_1,k_2)    if  k_1 ≤ m, s < i ≤ s+t
    
    C_i-s-t(k_1,k_2)    if  k_1 ≤ m, s+t < i ≤ s+2t
    
    D_i-s-2t(k_1,k_2)    if  k_1 ≤ m, s+2t < i ≤ r_CS
    Ŵ_i(k_1-m,k_2)    if  k_1 > m,  i ≤ r_CS
    Ẇ_i-r_CS(k_1,k_2)    if  r_CS < i 
    ,

for k_1=1,…, n, and k_2 = 1,…, m, where A_i := reshape(a_i,m × m ), for i=1, …, s, B_j:= reshape(b_j,m × m ), for j=1, …, t, and similarly for C_j and D_j, and 

    Ũ_k    := reshape( Ũ(:,k), m × (p-m) ), 
    Ṽ_k    := reshape( Ṽ(:,k),  p × (n-m) ),
    V̂_k    := reshape( V̂(:,k),  (p-m) × m ), 
    Ŵ_k    := reshape( Ŵ(:,k),  (n-m) × m ),

for k=1,…, r_CS, and

    U̇_i'   := reshape( U̇(:,i'), m × p ), 
    V̇_i'   := reshape( V̇(:,i'),  p × n ),
    Ẇ_i'   := reshape( Ẇ(:,i'),  n × m ),

for i'=1,…, (r-r_CS), and thus Ũ∈ℝ^m(p-m) × r_CS, Ṽ∈ℝ^p(n-m) × r_CS, V̂∈ℝ^(p-m)m × r_CS, Ŵ∈ℝ^(n-m)m × r_CS, U̇∈ℝ^mp × (r-r_CS), V̇∈ℝ^pn × (r-r_CS), and Ẇ∈ℝ^nm × (r-r_CS).
The proportions in <Ref> are those of T_234, (s,t):=(2,2), and r:=20. We show in the experiments that such PDs indeed exist. Remark that, contrary to the case when m, n, and p are equal, the CS rank r_CS does not equal the rank r, because this would be too restrictive. Also in the experiments that follow we did not find PDs for which both are equal. 
This is likely because of the interaction of the CS part with the parts Ũ, Ṽ, V̂, and Ŵ. However, including this structure still reduces the number of parameters from (mp+np+mn)r to (mp+np+mn)r-2m^2 r_CS.





§ NUMERICAL EXPERIMENTS
 

In this section we show the results that we obtained with the AL method from <cit.> to find PDs of MMTs with the generalized CS structure for different combinations of s and t. The AL method takes as input an upper and lower bound on the elements in the factor matrices: l ≤ x ≤ u.
To show the advantage of the new structure, we set u:=-l:=1 and generate 50 random starting points of size 10^-2 with the built-in function randn of :

    x_0:= 10^-2·randn((mp+np+mn)r-2m^2 r_CS× 1 ).

The number of inner iterations of the LM method was set to 50 and the number of outer iterations to 15. The tolerance on the gradient and on the bound constraint were both set to 10^-13.

In <Ref>, we show the results that we obtained with the AL method from <cit.> to find PD_11(T_223)s for different combinations of s and t. The second to last column (from left to right) indicates the number of exact solutions found for these 50 random starting points. With an exact solution, a solution for which the optimality conditions, i.e., the tolerance on the gradient and constraint, are met, and furthermore for which the cost function is smaller than 10^-12. The last column in the table indicates the number of practical solutions that were obtained when starting from the exact numerical PDs and by adding the constraint h_discr from (<ref>), scaled with a factor 0.1,
to the optimization problem. As can be seen, for many combinations, exact and practical PDs exist. An example practical PD_11(T_223) with (s,t):=(2,2) is given by the following equations:

    A    = [r,columns-width=0.4cm] 
    0     1 
    
    0     0 
    
    0     0 
    
    1     0
    , 
    B = [r,columns-width=0.4cm] 
     0     0
    
     0     0
    
    -1    -1
    
     0     0
    ,  C = [r,columns-width=0.4cm] 
    -1     0
    
     1    -1
    
    -1     0
    
     0     0
    ,  D = [r,columns-width=0.4cm]
     0     0
    
     1     1
    
     0    -1
    
     0     1
    ,
    V̂   = [r,columns-width=0.4cm]
    0     0     0     0     0     0     1    -1 
    
    0     0     0     0     0     0    -1     1
    , Ũ = [], Ṽ = [],  
    Ŵ   = [r,columns-width=0.4cm]
    0     0    -1     0     0     0     0    -1 
    
    0     0    -1     0     0     0     0    -1
    , 
    U̇   = [r,columns-width=0.4cm]
      0      1      0 
     
      0     -1     -1 
    
     -1      1      0 
    
      1     -1     -1 
    , Ẇ = [r,columns-width=0.4cm]
      0      0      0 
    
      0      0      1 
    
      0      0      0 
     
      0      0      0 
    
     -1     -1      0 
    
      1      0     -1 
    , V̇ = [r,columns-width=0.4cm]
      0      0      0 
    
      0      0      0 
    
      0     -1      1 
     
      1      0      0 
    
      0      0      0 
    
      1      0      1 
    .

Also the size and minimal, maximal, and number of different ranks of J_CS,gen, i.e., the Jacobian matrix of the new generalized CS structure, are given at the solutions. For a list of all the different ranks, we refer to the detailed results of the experiment, which are publicly available [<https://github.com/CharlotteVermeylen/ALM_FMM_gen_CS>]. Remark that the second dimension of the Jacobian matrix and the ranks decrease when r_CS increases. 



As can be seen from <Ref>, the highest value for r_CS that we were able to obtain is 10 and is obtained for (s,t):=(7,1). We were not able to transform this PD into a practical one using the constraint h_discr. 

In the last row, as a reference the results for (s,t):=(0,0), i.e., without the CS structure, are given. Remark that the number of numerical solutions in general decreases when r_CS increases but, on the other hand, the number of practical PDs increases for some combinations of s and t. 

In <Ref>, the convergence of the cost function is shown for the experiment with (s,t):=(7,0). As can be seen, 15 out of the 50 starting points converge quickly to a numerically exact PD in less than 100 iterations and the other starting points get stuck in swaps.



On the other hand, the results for PD_14(T_224)s are shown in <Ref>, again for different values of r_CS, s and t. The largest CS rank that we were able to obtain is again r_CS = 10 but now both PDs with (s,t)= (4,2) and (s,t)= (7,1) were found. However, again none of these PDs were discretizable using the constraint h_discr.

An example practical PD_14(T_224) that we obtained for (s,t)= (3,1) has as submatrices:

    A    = [r,columns-width=0.4cm] 
    1     1     0
    
    1     0     0
    
    0    -1     0
    
    0     0     1
    , 
    B = [r,columns-width=0.4cm] 
    -1 
    
    -1 
    
     1 
    
     1 
    ,  C = [r,columns-width=0.4cm] 
     0 
    
     0 
    
     1 
    
     0 
    ,  D = [r,columns-width=0.4cm] 
     0 
    
     1 
    
     0 
    
     0 
    ,
    Ṽ   = Ŵ = [r,columns-width=0.4cm]
    0     0     0     0     0     0 
    
    0     0     0     0     0     0 
    
    0     0     0     0     0     0 
    
    0     0     0     0     0     0 
    , Ũ = [], V̂ = [],
    U̇   = [r,columns-width=0.4cm]
     0     0     1     1     1     1     0    -1
    
    -1     0     0     1     1     0     0     0
    
     0    -1     0    -1     0     1     0     0
    
    -1    -1     0     0     0     0    -1     1
    ,
    V̇   = [r,columns-width=0.4cm]
    0     0     0    -1     0     0     0     0
    
    0     0     0    -1     0     0     0     0
    
    0     0     0     1     0     0     0     0
    
    0     0     0     0     0     0     0     0
    
    0     0    -1     0    -1     0     0     0
    
    0    -1     1     0     0    -1     0     1
    
    1     0     0     0     1     0     1    -1
    
    0     1     0     0     0     0    -1     0
    ,
    Ẇ   = [r,columns-width=0.4cm]
     0     0     0     1     0     0     0     0
    
     0     0     0     1     0     0     0     0
    
     0     0    -1     0     0    -1     0     0
    
     0    -1     0     0     0    -1    -1     1
    
     0     0     0    -1     0     0     0     0
    
     0     0     0     0     0     0     0     0
    
    -1     0     1     0    -1     0     0     1
    
    -1     0     0     0     0     0     1     0
    .

Again it can be noticed that by including the generalized CS structure, more practical PDs are obtained and furthermore, more different values of the rank of the Jacobian and thus more non-inv-equivalent PDs are obtained as long as the CS rank is not set too high.



Thirdly, in <Ref>, the results for T_233, rank 15 are shown. Now, the maximal value for r_CS that we obtained is 9, and one practical PD is obtained with parameters (s,t)=(6,1):

    A    = [r,columns-width=0.4cm] 
    0     1      0     0      0     0  
    
    0     0     -1     0      1     0  
    
    0     0      0     1     -1     0 
    
    0     0      1     1     -1     0
    , 
    B = [r,columns-width=0.4cm] 
     0 
    
     0 
    
     1 
    
     0 
    ,  C = [r,columns-width=0.4cm] 
      1 
    
      1 
    
     -1 
    
     -1
    ,  D = [r,columns-width=0.4cm] 
      0 
    
      1 
    
      0 
    
      0
    , 
    Ũ   = [r,columns-width=0.4cm] 
     1     0     0     0     0      1     0     0     0 
    
     1     0     0     0     0      0     0     0     0
    , 
    V̂   = [r,columns-width=0.4cm]
     0     0     0     0      0     1      0      0      1 
    
     0     1     0     1     -1    -1      1      1     -1 
    , 
    Ṽ   = [r,columns-width=0.4cm]
    -1     -1     0     0      0     0      0      0      1 
    
     0      0     1     1     -1     0     -1      0     -1 
    
    -1      0     0     0      0    -1      0      0      0
    , 
    Ŵ   = [r,columns-width=0.4cm]
     0     0     0     0     0     -1     0     0     0 
    
    -1     0     0     0     0     -1     0     0     0 
    ,
    U̇   = [r,columns-width=0.4cm]
      0      1      1      0      0     0 
    
     -1      0      1      1      1     1 
     
      0      1      0      0     -1     0 
    
      0      0      0      0      0     1 
    
      1     -1     -1      0      0     0 
     
      1      0     -1     -1     -1    -1 
    ,  
    V̇   = [r,columns-width=0.4cm]
      0      0      0      0      0      0 
    
      0      0      0      0      0      0 
    
     -1      0      0      1      0      0 
     
      0      0      0      0      0      0 
    
      0      0      0      0      0      0 
    
      1      1     -1     -1      1      0 
    
     -1      0      1      0      0      0 
     
      0      0      0      1     -1      1 
    
      0      0      0      0      0      0 
    , Ẇ = [r,columns-width=0.4cm]
     -1      0      1      1      1      0 
    
      0     -1      0      0      0      0 
    
     -1      1      1      1      1      0 
     
      0      0      0     -1     -1      0 
    
      0      0      0      0     -1     -1 
    
      1      0      0     -1      0      1 
    ,
 
 


Fourthly, in <Ref>, the results for T_225, rank 18, are shown. The maximal CS rank is 10 as in <Ref> and <Ref>. However, the highest value of the CS rank for which a practical PD is obtained is 9 and is obtained for (s,t) = (6,1). The CS submatrices of the PD that we obtained are:

    A    = [r,columns-width=0.4cm] 
    1    -1      1     0     0     0  
    
    0     1     -1     0     0     0  
    
    1    -1      0     0     0     0 
    
    0     0      0     1     1    -1
    , 
    B = [r,columns-width=0.4cm] 
     0 
    
     1 
    
     0 
    
     0 
    ,  C = [r,columns-width=0.4cm] 
      1 
    
     -1 
    
      1 
    
     -1
    ,  D = [r,columns-width=0.4cm] 
      0 
    
      0 
    
     -1 
    
      0
    ,


    Ṽ   = [r,columns-width=0.4cm]
     0     0     0     0     0     0     0     0     0 
    
     0     0     0    -1     0     1     0     0     0 
     
     0     0     0     0     0     0     0     0     0 
    
     0     0     0    -1     0     1     0     0     0 
    
     0     0     0     0     0     0     0     0     0 
     
     0     0     0     0    -1     1     0     0     0 
    , Ũ = [], V̂ = [], 
    
    Ŵ   = [r,columns-width=0.4cm]
     0     0     0     0     0     0     0     0     0 
    
     0     0     0     0    -1     1     0     0     0 
     
     0     0     0    -1     0     1     0     0     0 
    
     0     0     0     0     0     0     0     0     0 
    
     0     0     0     0     1    -1     0     0     0 
     
     0     0     0     1     0    -1     0     0     0 
    ,
    
    U̇   = [r,columns-width=0.4cm]
     0     1     0     1     0     0    -1    -1     1 
    
     0     0     0     0     1     1     0     0     1 
     
     1     0     1    -1     0     0     1     0    -1 
    
     1     0     1    -1     0     0     0     0    -1 
    , 
    
    V̇   = [r,columns-width=0.4cm]
     0     0     0     0     0     0     0     0     0 
    
     0     0     0     0     0     0     0     0     0 
     
     0     0     0     0     0     0     0     0     0 
    
     0     0     0     0     0     0     0     0     0 
    
     0     1     0     0    -1     0     0     0     0 
     
     0     1     0    -1     0     0    -1     0     0 
    
     0     0     1     1    -1     0     0     0     1 
    
     0     0     1     0     0     0     0     0     0 
     
     0     0     1     1     0    -1     0    -1     1 
    
    -1     0     1     0     0     0     0     0     0 
    ,
    
    Ẇ   = [r,columns-width=0.4cm]
     0     0     0     0     0     0     0     0     0 
    
     0     0     0     0     0     0     0     0     0 
     
     0     1     0     0     0     0    -1     0     0 
    
     1     0     1     1     0     0     1    -1     0 
    
    -1     0     0     0     0     0     0     1     0 
     
     0     0     0     0     0     0     0     0     0 
    
     0     0     0     0     0     0     0     0     0 
    
     0     0     0     1    -1    -1     1     0    -1 
     
     0     0     0    -1     0     1    -1     0     1 
    
     0     0     0     0     0    -1     0     0     0 
    .





























































Then in <Ref>, the results for T_234, rank 20, are shown. The parameters s and t for which solutions are found are similar to the previous experiments but the number of solutions is on average smaller because the problem becomes more difficult as p, n, and the rank r increase.
Also, fewer different values of the rank of the Jacobian matrix are found. Remark that the number of different values of the rank in general increases for moderate values of r_CS. Furthermore, including the CS structure again helps to find practical solutions compared to the generic case ((s,t)=(0,0)) for different combinations of s and t. 



Lastly in <Ref>, the results for T_245 and rank 33 are shown. Again the problem has become more difficult because p, n, and r have increased. 
Still, we find more numerical solutions for several CS parameters compared to the generic case ((s,t) = (0,0)). 





Remark that we can also choose the CS rank smaller than r when m=p=n. For example, we are able to enforce 11 symmetric rank-1 tensors in a PD_23(T_333). The other 12 rank-1 tensors do not have to admit a structure in this case. One of the practical PDs that we found for these parameters has as symmetric part:

    A    = [r,columns-width=0.4cm]
     1     0     0     0     0     0     0     0     0     0     0
    
     0     0     0     1     0     0     0     0    -1     0     0
    
     0    -1    -1     0     1     0     0     0     0     0     1
    
     0     0    -1     0     0     1    -1    -1     1     0     1
    
     0     0     0     1     0     0     1     0    -1     0     0
    
     0    -1    -1     0     1     1     0    -1     0     0     1
    
     0     0     0     0     0     0     0     1     0     0    -1
    
     0     0     1     0    -1    -1     0     0     0     1     0
    
     0     1     1     0    -1    -1     0     1     0     1    -1
    ,

and as asymmetric part:

    U̇   = [r,columns-width=0.4cm]
    -1     0     0     0     0    -1     0     0     0     0     0     0
    
     1     0     0     1     0     0    -1     0    -1     0     0     0
    
     0     0     0     1     0    -1     0     0    -1    -1     0     0
    
    -1     0     0     0     0    -1     0     0     0     0    -1     1
    
     1     0     0     1     0     0     0     0     0     0     0     0
    
     0     0     0     1     0    -1     0     1     0     0    -1     0
    
     1     0     1     0     1     1     0     0     0     0     0     0
    
    -1     1    -1    -1     0     0     0     0     0     0     0     0
    
     0     0     0    -1     0     1     0     0     0     0     0     0
    , 
    V̇   = [r,columns-width=0.4cm]
     0     0     0     0     0     0    -1     0     0     1     0     0
    
     0     0     0     0     0     0     1     1     0     0    -1     1
    
     0     0     0     0     1     0     0     1     0     1    -1     0
    
     1     0    -1     0     0     0    -1     0     0     1     0     0
    
     0     0     0     0     0     0     1     1     0     0     0     0
    
     0     1    -1     0     0     0     0     1     0     1     0     0
    
     0     0     0     0     0    -1     1     0     1    -1     0     0
    
     0     0     0    -1     0     0    -1    -1    -1     0     0     0
    
     0     0     0     0     0     0     0    -1     0    -1     0     0
    ,


    Ẇ   = [r,columns-width=0.4cm]
     0     0     0     0     1     0     0     0     0     0     0     1
    
    -1     1    -1     0     0     0     0     0     0     0     0    -1
    
     0     1    -1     0     1     1     0     0     0     0     0     0
    
     0     0     0     0     1     0     1     0    -1     0     0     1
    
     0     1     0     0     0     0     0     0     0     0     0    -1
    
     0     1     0    -1     1     0     0     0    -1     0     0     0
    
     0     0     0     0    -1     0     0     0     0    -1     1    -1
    
     0    -1     0     0     0     0     0     1     0     0    -1     1
    
     0    -1     0     0    -1     0     0     0     0     0     0     0
    .






§ CONCLUSION

The proposed generalized cyclic symmetric structure can be useful to include in the optimization problem for the fast matrix multiplication (FMM) problem to reduce the number of parameters and to improve the convergence to numerically exact and practical polyadic decompositions (PDs) of matrix multiplication tensors (MMTs). Several new practical PDs, i.e., sparse PDs with elements that are either one, minus one, or zero, with this new structure are discovered for various problem parameters, i.e., sizes of the MMT. 



§ ACKNOWLEDGEMENT


This work was funded by (1) the Flemish Government: this research received funding under the AI Research Program. Charlotte Vermeylen is affiliated to Leuven.AI - KU Leuven institute for AI, B-3000, Leuven, Belgium. (2) KU Leuven Internal Funds: iBOF/23/064, C14/22/096 and
IDN/19/014.
Marc Van Barel was supported by the Research Council KU Leuven, C1-project C14/17/073 and by the Fund for Scientific Research–Flanders (Belgium), EOS Project no 30468160 and project G0B0123N.



§ COMPETING INTERESTS


Declarations of interest: none.










elsarticle-num 














