Equivalent Polyadic Decompositions of Matrix Multiplication Tensors

TL;DR

This paper introduces an efficient algorithm to determine the equivalence of polyadic decompositions of matrix multiplication tensors, analyzes their classes for various sizes, and explores conditions for integer-based decompositions to improve fast matrix multiplication algorithms.

Contribution

It presents a novel algorithm for deciding equivalence of tensor decompositions and analyzes their classes for different matrix sizes, advancing understanding of fast matrix multiplication methods.

Findings

For larger matrices, decompositions are likely essentially different.

The algorithm effectively distinguishes equivalent decompositions.

Necessary conditions for integer-based decompositions are identified.

Abstract

Invariance transformations of polyadic decompositions of matrix multiplication tensors define an equivalence relation on the set of such decompositions. In this paper, we present an algorithm to efficiently decide whether two polyadic decompositions of a given matrix multiplication tensor are equivalent. With this algorithm, we analyze the equivalence classes of decompositions of several matrix multiplication tensors. This analysis is relevant for the study of fast matrix multiplication as it relates to the question of how many essentially different fast matrix multiplication algorithms there exist. This question has been first studied by de~Groote, who showed that for the multiplication of $2\times2$ matrices with $7$ active multiplications, all algorithms are essentially equivalent to Strassen's algorithm. In contrast, the results of our analysis show that for the multiplication of larger matrices, (e.g., $2\times3$ by $3\times2$ or $3\times3$ by $3\times3$ matrices), two decompositions are very likely to be essentially different. We further provide a necessary criterion for a polyadic decomposition to be equivalent to a polyadic decomposition with integer entries. Decompositions with specific integer entries, e.g., powers of two, provide fast matrix multiplication algorithms with better efficiency and stability properties. This condition can be tested algorithmically and we present the conclusions obtained for the decompositions of small/medium matrix multiplication tensors.

Figures

Tables

Table 1. Table 1: First and second columns: different cases considered in the numerical experiments. Fourth column: total time required to compute the N s = 10 000 subscript 𝑁 𝑠 10000 N_{s}=10\,000 decompositions with Tichavský et al.’s method [ 26 ] . Third column: number of randomly generated initial guesses (trials) ( 𝐔 [ F ] , 0 , 𝐕 [ F ] , 0 , 𝐖 [ F ] , 0 ) subscript 𝐔 delimited-[] 𝐹 0 subscript 𝐕 delimited-[] 𝐹 0 subscript 𝐖 delimited-[] 𝐹 0 (\mathrm{\mathbf{U}}_{{[F]},0},\mathrm{\mathbf{V}}_{{[F]},0},\mathrm{\mathbf{W}}_{{[F]},0}) we had to use to compute the N s = 10 000 subscript 𝑁 𝑠 10000 N_{s}=10\,000 samples.

$(m,p,n)$	$F$	# trials	Elapsed time [hours]
$(1,2,1)$	$2$	$\mathrm{10\hspace{0.17em}000}$	$0.017$
$(2,1,2)$	$4$	$\mathrm{10\hspace{0.17em}000}$	$0.05$
$(2,2,2)$	$7$	$\mathrm{10\hspace{0.17em}762}$	$0.32$
$(2,3,2)$	$11$	$\mathrm{10\hspace{0.17em}133}$	$0.51$
$(3,2,3)$	$15$	$\mathrm{18\hspace{0.17em}003}$	$40$
$(3,3,3)$	$23$	$\mathrm{15\hspace{0.17em}829}$	$16.8$

Table 2. Table 2: Equivalence of decompositions. For each case, we have used Algorithm  1 to check the equivalence between 10 000 10000 10\,000 randomly chosen pairs of decompositions inside the cluster. The second column gives the percentage of pairs of equivalent decompositions. The third column gives the average time required to check the equivalence of the decompositions with Algorithm  1 . ∗ The depth of the algorithm is the maximal length of a partial permutation π ∈ Inj ​ ( n , F ) 𝜋 Inj 𝑛 𝐹 \pi\in\mbox{Inj}(n,F) that is rejected (see Algorithm  1 ).

	Percentage of equivalent pairs	Mean elapsed time [sec]	Depth∗
			Max.	Mean
$(1,2,1)$	$100\%$	$4.74\cdot {10}^{-4}$	1	1
$(2,1,2)$	$0\%$	$1.02\cdot {10}^{-3}$	1	1
$(2,2,2)$	$100\%$	$2.68\cdot {10}^{-3}$	4	0.55
$(2,3,2)$	$0\%$	$4.94\cdot {10}^{-3}$	6	1.08
$(3,2,3)$	$0\%$	$2.66\cdot {10}^{-2}$	9	2.91
$(3,3,3)$	$0\%$	$2.82\cdot {10}^{-2}$	5	1.51

Equations

$${\Phi }_{m,p,n}(\mathbf{A},\mathbf{B})=\mathbf{AB}.$$

$$\Phi (u,v)=\sum _{r=1}^F{f}_r(u)g_r(v)w_r$$

$$f_r(\mathbf{A})=\mathrm{trace}({\mathbf{U}}_r\mathbf{A})=\sum _{i=1}^p\sum _{j=1}^m{\mathbf{U}}_r^{(i,j)}{\mathbf{A}}^{(j,i)}$$

$${\mathbf{U}}_r^{\prime }={\mathbf{U}}_{\sigma (r)},{\mathbf{V}}_r^{\prime }={\mathbf{V}}_{\sigma (r)},{\mathbf{W}}_r^{\prime }={\mathbf{W}}_{\sigma (r)}.$$

$${\mathbf{U}}_r^{\prime }={\lambda }_r{\mathbf{U}}_r,{\mathbf{V}}_r^{\prime }={\mu }_r{\mathbf{V}}_r,{\mathbf{W}}_r^{\prime }={\nu }_r{\mathbf{W}}_r.$$

$${\mathbf{U}}_r^{\prime }={\mathbf{Q}}^{-1}{\mathbf{U}}_r\mathbf{P},{\mathbf{V}}_r^{\prime }={\mathbf{R}}^{-1}{\mathbf{V}}_r\mathbf{Q},{\mathbf{W}}_r^{\prime }={\mathbf{P}}^{-1}{\mathbf{W}}_r\mathbf{R}.$$

$$f_r^{\prime }(\mathbf{A})=\mathrm{trace}({\mathbf{U}}_r^{\prime }\mathbf{A})=\mathrm{trace}({\mathbf{U}}_r[\mathbf{PA}{\mathbf{Q}}^{-1}])=f_r(\mathbf{PA}{\mathbf{Q}}^{-1})$$

$${\Phi }^{\prime }(\mathbf{A},\mathbf{B})$$

$$={\mathbf{P}}^{-1}{\Phi }_{m,p,n}(\mathbf{PA}{\mathbf{Q}}^{-1},\mathbf{QB}{\mathbf{R}}^{-1})\mathbf{R}={\mathbf{P}}^{-1}(\mathbf{PA}{\mathbf{Q}}^{-1})(\mathbf{QB}{\mathbf{R}}^{-1})\mathbf{R}=\mathbf{AB}.$$

$$\mathbf{P}={\left(\frac{\lambda }{\nu }\right)}^{1/3}{\mathbf{I}}_m,\mathbf{Q}={\left(\frac{\mu }{\lambda }\right)}^{1/3}{\mathbf{I}}_p,\mathbf{R}={\left(\frac{\nu }{\mu }\right)}^{1/3}{\mathbf{I}}_n$$

$$\mathbf{MA}=\mathbf{A}\mathrm{diag}({\xi }_1,\dots ,{\xi }_n)$$

$$\mathbf{M}={\mathbf{A}}^{\prime }\mathrm{diag}({\xi }_1,\dots ,{\xi }_m){\mathbf{A}}^{\prime -1}$$

$$\mathrm{\mathbf{M}}\mathrm{\mathbf{a}}_{j}=\big{[}\mathrm{\mathbf{A}}^{\prime}\,\mathrm{diag}(\xi_{1},\ldots,\xi_{m})\,\mathrm{\mathbf{A}}^{\prime-1}\big{]}\mathrm{\mathbf{a}}_{j}=\sum_{i=1}^{m}\mathrm{\mathbf{a}}_{i}\xi_{i}q_{j,i}=\mathrm{\mathbf{a}}_{j}\xi_{j}=\xi_{j}\sum_{i=1}^{m}\mathrm{\mathbf{a}}_{i}q_{j,i}.$$

$${\mathbf{M}}_s(\mathbf{x}+\mathbf{y})=\mathbf{x}\text{for all x∈Us and y∈⨁s′=sUs′}.$$

$$\dim (\mathcal{S})={\mathrm{cl}}_{\oplus }(\mathbf{A})+(m-1)(m-r)+Z.$$

$$\mathbf{MXA}=\mathbf{XA}\mathrm{diag}({\xi }_1,\dots ,{\xi }_n)$$

$$\mathrm{\mathbf{M}}=\leavevmode\hbox to83.51pt{\vbox to83.1pt{\pgfpicture\makeatletter\hbox{\hskip 240.92488pt\lower-84.23367pt\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ } {}{}{}{{}}{{}}{{}}{ {}}{{}}\hbox{\hbox{\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{\offinterlineskip{}{}{{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}}{{{}}}{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-233.16931pt}{-73.83368pt}\pgfsys@invoke{ }\hbox{\vbox{\halign{\pgf@matrix@init@row\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding&&\pgf@matrix@step@column{\pgf@matrix@startcell#\pgf@matrix@endcell}&#\pgf@matrix@padding\cr\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip 10.0pt\hfil&\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip 10.0pt\hfil&\hfil\hskip 4.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{}}}&\hskip 4.0pt\hfil&\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.0pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.0pt\hfil&\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.0pt\hfil&\hfil\hskip 4.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 4.0pt\hfil&\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.0pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip 4.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 4.0pt\hfil&\hfil\hskip 4.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 4.0pt\hfil&\hfil\hskip 4.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 4.0pt\hfil&\hfil\hskip 4.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 4.0pt\hfil\cr\vskip 0.0pt\cr\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.0pt\hfil&\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.0pt\hfil&\hfil\hskip 4.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{$$

} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 4.0pt\hfil&\hfil\hskip 10.0pt\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}}\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{0.0pt}{0.0pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{

$$} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}&\hskip 10.0pt\hfil\cr}}}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}{{{{}}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ }}{ } {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-240.92488pt}{-73.53368pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{{\nullfont{{{ {}{}{}}}{{ {}{}{}}}}}$\left[\vbox{\hrule height=36.2pt,depth=36.2pt,width=0.0pt}\right.$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}}\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{}}{} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-162.96931pt}{-73.53368pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{{\nullfont{{{ {}{}{}}}{{ {}{}{}}}}}$\left.\vbox{\hrule height=36.2pt,depth=36.2pt,width=0.0pt}\right]$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{{}}{} { {}}{}{{}}{}{}{}{}{{}}{}{}{}{{{}{}}}{}\pgfsys@moveto{-233.36931pt}{-15.63368pt}\pgfsys@moveto{-233.36931pt}{-15.63368pt}\pgfsys@lineto{-233.36931pt}{-56.03368pt}\pgfsys@lineto{-192.96931pt}{-56.03368pt}\pgfsys@lineto{-192.96931pt}{-15.63368pt}\pgfsys@closepath\pgfsys@moveto{-192.96931pt}{-56.03368pt}\pgfsys@stroke\pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-219.15265pt}{-38.34811pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\mathrm{\mathbf{M}}_{1}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{{}}{} { {}}{}{{}}{}{}{}{}{{}}{}{}{}{{{}{}}}{}\pgfsys@moveto{-233.36931pt}{-63.63368pt}\pgfsys@moveto{-233.36931pt}{-63.63368pt}\pgfsys@lineto{-233.36931pt}{-84.03368pt}\pgfsys@lineto{-192.96931pt}{-84.03368pt}\pgfsys@lineto{-192.96931pt}{-63.63368pt}\pgfsys@closepath\pgfsys@moveto{-192.96931pt}{-84.03368pt}\pgfsys@stroke\pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-219.15265pt}{-76.34811pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\mathrm{\mathbf{M}}_{2}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} { {}}{}{{}}{} { {}}{}{{}}{}{}{}{}{{}}{}{}{}{{{}{}}}{}\pgfsys@moveto{-185.36931pt}{-15.63368pt}\pgfsys@moveto{-185.36931pt}{-15.63368pt}\pgfsys@lineto{-185.36931pt}{-84.03368pt}\pgfsys@lineto{-164.96931pt}{-84.03368pt}\pgfsys@lineto{-164.96931pt}{-15.63368pt}\pgfsys@closepath\pgfsys@moveto{-164.96931pt}{-84.03368pt}\pgfsys@stroke\pgfsys@invoke{ }\hbox{\hbox{{\pgfsys@beginscope\pgfsys@invoke{ }{{}{}{{ {}{}}}{ {}{}} {{}{{}}}{{}{}}{}{{}{}} { }{{{{}}\pgfsys@beginscope\pgfsys@invoke{ }\pgfsys@transformcm{1.0}{0.0}{0.0}{1.0}{-181.15265pt}{-52.34811pt}\pgfsys@invoke{ }\hbox{{\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{ }\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{ }\hbox{{$\mathrm{\mathbf{M}}_{3}$}} }}\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope{{ {}{}{}}{{{}}{{}}}}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\hskip 28.45274pt\text{with}\;\;\;\left\{\begin{array}[]{r@{\:}c@{\:}l}\mathrm{\mathbf{M}}_{1}\>&\in\hfil\>&\mathbb{R}^{r\times r}\\ \mathrm{\mathbf{M}}_{2}\>&\in\hfil\>&\mathbb{R}^{(m-r)\times r}\\ \mathrm{\mathbf{M}}_{3}\>&\in\hfil\>&\mathbb{R}^{m\times(m-r)}\\ \end{array}\right.\;.$$

$${\mathbf{M}}_1\tilde{\mathbf{A}}=\tilde{\mathbf{A}}\mathrm{diag}({\xi }_1,\dots ,{\xi }_n)\text{and}{\mathbf{M}}_2\tilde{\mathbf{A}}=0.$$

$$\mathrm{trace}(\mathbf{B}{\mathbf{Z}}^{\top })$$

$$=\sum _{r=1}^F\mathrm{trace}({\mathbf{U}}_r{\mathbf{A}}_{\ast })\mathrm{trace}({\mathbf{V}}_r\mathbf{B})\mathrm{trace}({\mathbf{W}}_r\mathbf{C}).$$

$$\tilde{\mathbf{U}}=[\mathrm{vec}({\mathbf{U}}_1),\dots ,\mathrm{vec}({\mathbf{U}}_F)]\in {\mathbb{R}}^{pm\times F}$$

$${\lambda }_r{\mathbf{U}}_r^{\prime }={\mathbf{Q}}^{-1}{\mathbf{U}}_r\mathbf{P},{\mu }_r{\mathbf{V}}_r^{\prime }={\mathbf{R}}^{-1}{\mathbf{V}}_r\mathbf{Q},{\nu }_r{\mathbf{W}}_r^{\prime }={\mathbf{P}}^{-1}{\mathbf{W}}_r\mathbf{R}$$

$${\Phi }^{\prime \prime }(u,v)=\sum _{r=1}^F{\lambda }_r{\mu }_r{\nu }_r{f}_r^{\prime }(u)g_r^{\prime }(v)w_r^{\prime }.$$

$$\mathrm{vec}(\mathbf{AXB})=({\mathbf{B}}^{\top }\otimes \mathbf{A})\mathrm{vec}(\mathbf{X}).$$

$$({\mathbf{P}}^{\top }\otimes {\mathbf{Q}}^{-1})\tilde{\mathbf{U}}={\tilde{\mathbf{U}}}^{\prime }\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_F).$$

$$\mathbf{M}\tilde{\mathbf{U}}={\tilde{\mathbf{U}}}^{\prime }\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_F)$$

$$\mathbf{MA}={\mathbf{A}}^{\prime }\mathrm{diag}({\xi }_1,\dots ,{\xi }_n)$$

$$\left\{\begin{array}[]{r@{\;\,}c@{\;\,}ll}\mathrm{\mathbf{R}}\mathrm{\mathbf{V}}^{\prime}_{r}&=&\tilde{\mu}_{r}\mathrm{\mathbf{V}}_{r}\mathrm{\mathbf{Q}}&\quad\text{for all $1\leq r\leq F$},\\[2.0pt] \mathrm{\mathbf{W}}_{r}^{\prime}\mathrm{\mathbf{R}}&=&\tilde{\nu}_{r}\mathrm{\mathbf{P}}\mathrm{\mathbf{W}}_{r}&\quad\text{for all $1\leq r\leq F$},\end{array}\right.$$

$$\begin{array}{rl}{\mathbf{M}}_r^{\prime }& ={\lambda }_r{\mu }_r{\nu }_r{\mathbf{M}}_r^{\prime }\\ & =({\mathbf{P}}^{-1}{\mathbf{W}}_{\sigma (r)}\mathbf{R})({\mathbf{R}}^{-1}{\mathbf{V}}_{\sigma (r)}\mathbf{Q})({\mathbf{Q}}^{-1}{\mathbf{U}}_{\sigma (r)}\mathbf{P})={\mathbf{P}}^{-1}{\mathbf{M}}_{\sigma (r)}\mathbf{P}.\end{array}$$

$$\pi\big{(}(\mathrm{\mathbf{M}}_{r})_{r\in I}\big{)}\;\;\;\text{and}\;\;\;(\mathrm{\mathbf{M}}^{\prime}_{r})_{r\in I}\;\;\;\text{are simultaneously similar.}$$

$$\mathbf{X}{\mathbf{A}}_i-{\mathbf{B}}_i\mathbf{X}=0\text{for all}1\le i\le m,$$

Full text

††journal: Journal of Computational and Applied Mathematics

1 Introduction

The straightforward way to multiply two $N\times N$ matrices costs $\mathcal{O}(N^{3})$ operations. In particular, multiplying two $2\times 2$ matrices requires $8$ scalar multiplications. However, as first remarked by V. Strassen [1] in 1969, the arithmetic operations can be grouped cleverly to reduce the work to $7$ multiplications only. By doing this recursively, we can reduce the cost for the multiplication of $N\times N$ matrices to $\mathcal{O}(N^{2.81})$ operations. Strassen’s discovery opened the door to a considerable amount of research on the algorithmic complexity of matrix multiplication (see paragraphs below for a survey). The reduction of the complexity may actually become so significant that a new architecture for large matrix multiplication is emerging. Essential is first that we find inexpensive schemes for the multiplication of relatively small matrices.

The multiplication of $m\times p$ matrices by $p\times n$ matrices can be represented by a third-order tensor. Finding inexpensive schemes for the multiplication of such matrices can be approached by decomposing the associated tensor as a sum of rank-$1$ terms (such a decomposition is called polyadic decomposition). The minimal number of rank-$1$ terms necessary to decompose a tensor is its rank. In the case of matrix multiplication, the rank of the associated tensor is equal to the smallest number of active multiplications needed to compute the matrix product. (By active multiplication, we mean a multiplication of two scalars that both depend on the matrices to be multiplied.) As a consequence, determining the rank of the associated tensor allows us to find an exponent $\alpha$ such that the complexity for the multiplication of $N\times N$ matrices is at most of $\mathcal{O}(N^{\alpha})$ arithmetic operations [2].

Although the problem of matrix multiplication complexity is quite old, only partial results are known so far. Even for the multiplication of small matrices, determining the rank of the associated tensor is still an open problem. The largest case that is completely understood is the multiplication of $2\times 2$ matrices by $2\times 2$ matrices. The rank of the associated tensor is $7$ [3] (so that Strassen’s algorithm is optimal), and it was proved by de Groote [4] that the decomposition induced by Strassen’s algorithm is essentially unique (with respect to a class of transformations acting on polyadic decompositions of matrix multiplication; see the paragraph hereunder). For the multiplication of $3\times 3$ matrices, an algorithm using $23$ active multiplications was proposed by Laderman [5] in 1976, and by Makarov in 1986 [6]; see also [7]. This means that the rank of the associated tensor is at most $23$. On the other hand, Bläser proved [8] in 2003 that the rank for the multiplication of $3\times 3$ matrices should be at least $19$. The gap $19$–$23$ has not been reduced since then. Other algorithms for the multiplication of matrices with small and medium sizes are proposed in [9], [10], [11], [12]; also the paper [13] of Ballard et al. gives a good overview of the practical algorithms that are available in 2016. Further reductions of the exponent of matrix multiplication complexity have been achieved, e.g., by Pan [14, 15], Bini et al. [16], Schönhage [17], and Coppersmith and Winograd [18] by means of more advanced techniques (including namely the study of the “border rank” of the associated tensor; see also [19, 20]). Currently, the best known upper bound for the complexity of matrix multiplication is $\mathcal{O}(N^{2.3729})$ by Le Gall [21].

The aim of this paper is to study the connections between polyadic decompositions for matrix multiplication tensors. In particular, we consider three types of transformations, called invariance transformations, acting on the set of polyadic decompositions of a given matrix multiplication tensor, and we study the equivalence relation induced by these transformations. These transformations have been studied by de Groote [22, 4] who has shown that Strassen’s algorithm is essentially unique in the sense that every other decomposition with $7$ rank-$1$ terms is equivalent to it. In contrast, for the multiplication of $3\times 3$ matrices, Johnson and McLoughlin [23] showed that Laderman’s algorithm [5] is not essentially unique. They provided two parametrized families of decompositions (with $23$ rank-$1$ terms) of the $3\times 3$ matrix multiplication tensor that are mutually inequivalent and also inequivalent to Laderman’s. Later, Oh, Kim and Moon [24] discovered other decompositions of the $3\times 3$ matrix multiplication tensor inequivalent to the previous ones, and Sedoglavic [7] showed that Laderman’s algorithm can be constructed using Strassen’s algorithm and related tensor’s transformations.

The techniques used by de Groote, Johnson and McLoughlin, and Oh et al. to prove the equivalence/inequivalence of decompositions are either too specific [4, 23] or too conservative [24] (some inequivalent decompositions are not recognized as such) to be applied to general decompositions of matrix multiplication tensors of arbitrary size. In this paper, we present an algorithm for deciding whether two given decompositions are equivalent through invariance transformations. Thanks to this algorithm, we were able to study the equivalence classes of large sample sets of matrix multiplication tensor decompositions (computed with numerical methods, see Section 6). This allows us to get a better understanding of the equivalence relation of decompositions: for instance, the numerical experiments (Section 6.3) suggest that for tensors larger than the $2\times 2$ by $2\times 2$ case, two “generic” decompositions are inequivalent.

In addition, we describe a necessary criterion for a matrix multiplication tensor decomposition to be discretizable, that is, to be equivalent to a decompositions whose rank-$1$ terms can be factorized into vectors or matrices whose entries only take a few distinct values (for instance, we may want that all entries of the factor vectors/matrices of the rank-$1$ terms belong to the set $\{-1,0,+1\}$). Such decompositions are called discrete decompositions [25]. Our interest in discretizable decompositions originates from the observation that for small/medium matrix multiplication, the decompositions proposed in the literature are generally discrete: Strassen’s and Laderman’s algorithms are discrete with factor matrices coefficients belonging to $\{-1,0,+1\}$, other (inequivalent) decompositions with coefficients in $\{-1,0,+1\}$ are proposed, e.g., in [13, 24, 26] for the multiplication of $2\times 3$ by $3\times 2$, $2\times 3$ by $3\times 3$, and $3\times 3$ by $3\times 3$ matrices. In particular, all the decompositions listed in [13] are discrete.

Discrete decompositions provide matrix multiplication algorithms with better efficiency and stability properties. However, the classical iterative processes for computing tensor decompositions do not lead in general to solutions of this kind. A reasonable approach to compute discrete decompositions is then to (i) compute a general decomposition, and (ii) use invariance transformations to obtain an equivalent discrete decomposition. Closely related methods are used, e.g., in [23, 26]. The necessary criterion for discretizability allows us to identify some decompositions that cannot be transformed via invariance transformations into a discrete decomposition with a given “target set” for the coefficients. By applying the necessary criterion to the sample sets of decompositions, we observed that, contrary to what the decompositions available in the literature suggest, most of the decompositions for tensors larger than the $2\times 2$ by $2\times 2$ case are not discretizable with respect to the commonly-used target sets (e.g., $\{0,\pm 1\}$ or $\{0,\pm 1/2,\pm 1\}$).

The paper is organized as follows. In Section 2, we introduce the notation and recall the definitions of matrix multiplication tensors and polyadic decompositions. Invariance transformations and the induced equivalence relations are introduced in Section 3. In Section 4, we describe the algorithm for deciding whether two decompositions are equivalent and if so, computing the invariance transformations involved in their equivalence. The necessary criterion for discretizability is discussed in Section 5. Numerical experiments are presented in Section 6.

2 Preliminaries

2.1 Matrix multiplication tensors and polyadic decompositions

Let $U$, $V$ and $W$ be vector spaces over a field $\mathbb{F}$. We denote by $\mathrm{Bil}(U,V;W)$ the set of $\mathbb{F}$-bilinear maps from $U\times V$ to $W$. For positive integers $m,p,n$, the multiplication of $m\times p$ matrices by $p\times n$ matrices can be represented by the bilinear map $\Phi_{m,p,n}\in\mathrm{Bil}(\mathbb{F}^{m\times p},\mathbb{F}^{p\times n};\mathbb{F}^{m\times n})$ defined by

[TABLE]

From the identification between multilinear maps and tensors, $\Phi_{m,p,n}$ is sometimes referred to as the $(m,p,n)$ matrix multiplication tensor.

The concept of rank of a bilinear map $\Phi\in\mathrm{Bil}(U,V;W)$ is central in the analysis of the asymptotic complexity of matrix multiplication. We say that $\Phi\neq 0$ has rank $1$ if $\Phi(u,v)=f(u)g(v)w$ for some $f\in U^{*}$, $g\in V^{*}$ and $w\in W$, where $U^{*}$ and $V^{*}$ are the dual spaces of $U$ and $V$ respectively. For a general $\Phi\in\mathrm{Bil}(U,V;W)$, an $F$-term polyadic decomposition (in short $F$-PD) of $\Phi$ is a decomposition of $\Phi$ as the sum of $F$ rank-$1$ terms [27, 28]:

[TABLE]

for some $f_{r}\in U^{*}$, $g_{r}\in V^{*}$ and $w_{r}\in W$. The rank of $\Phi$ is the smallest $F$ such that $\Phi$ admits an $F$-term polyadic decomposition (1).

For a matrix multiplication tensor $\Phi_{m,p,n}$, a polyadic decomposition like (1) requires $f_{r}\in(\mathbb{F}^{m\times p})^{*}$, $g_{r}\in(\mathbb{F}^{p\times n})^{*}$ and $w_{r}=\mathrm{\mathbf{W}}_{r}\in\mathbb{F}^{m\times n}$. We may identify $f_{r}$ with the unique matrix $\mathrm{\mathbf{U}}_{r}\in\mathbb{F}^{p\times m}$ such that

[TABLE]

for every $\mathrm{\mathbf{A}}\in\mathbb{R}^{m\times p}$, where $\mathrm{\mathbf{M}}^{(i,j)}$ denotes the $(i,j)$th entry of a matrix $\mathrm{\mathbf{M}}$. In the same way, $g_{r}$ can be identified with a unique matrix $\mathrm{\mathbf{V}}_{r}\in\mathbb{F}^{n\times p}$. If $\mathrm{\mathbf{U}}_{1},\ldots,\mathrm{\mathbf{U}}_{F}$, $\mathrm{\mathbf{V}}_{1},\ldots,\mathrm{\mathbf{V}}_{F}$ and $\mathrm{\mathbf{W}}_{1},\ldots,\mathrm{\mathbf{W}}_{F}$ give rise to a decomposition (1) of $\Phi_{m,p,n}$, we will say with slight abuse of notation that the triple $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$444The rational behind this notation is that if $\mathrm{\mathbf{X}}_{1},\mathrm{\mathbf{X}}_{2},\ldots,\mathrm{\mathbf{X}}_{F}$ are $F$ mathematical objects, then $\mathrm{\mathbf{X}}_{[F]}$ denotes the ordered set, or $F$-uple, $(\mathrm{\mathbf{X}}_{1},\mathrm{\mathbf{X}}_{2},\ldots,\mathrm{\mathbf{X}}_{F})$. is an $F$-term polyadic decomposition ($F$-PD) of $\Phi_{m,p,n}$.

The link between matrix multiplication complexity and the rank of matrix multiplication tensors is nicely explained in [2, Chapter 15]. Especially, it is shown how to build from an $F$-term polyadic decomposition ($F$-PD) of $\Phi_{m,p,n}$ a recursive algorithm for the multiplication of $N\times N$ matrices over $\mathbb{F}$ with complexity in $\mathcal{O}(N^{\omega+\varepsilon})$ arithmetic operations $\{+,-,\times\}$, with $\omega=3\log_{mpn}(F)$ and for any $\varepsilon>0$. For instance, Strassen’s algorithm can be obtained from a decomposition of $\Phi_{2,2,2}$ with $7$ terms. This directly gives the well-known upper bound $\omega=3\log_{8}(7)\approx 2.81$ for the exponent of matrix multiplication complexity [1].

2.2 Discrete decompositions

In this paper, we focus on algorithms for matrix multiplication over the field of real numbers, i.e., on the case $\mathbb{F}=\mathbb{R}$ and $\mathrm{\mathbf{U}}_{r}$, $\mathrm{\mathbf{V}}_{r}$ and $\mathrm{\mathbf{W}}_{r}$ are real matrices.

For the problem of matrix multiplication over $\mathbb{R}$ (or $\mathbb{C}$), two decompositions might be not equally useful even if they have the same number of rank-$1$ terms. For instance, decompositions with “structured” values in the rank-$1$ terms are more useful in practice. This leads us to the following definition:

Definition 2.1**.**

A decomposition $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ of $\Phi_{m,p,n}$ is said to be discrete if the entries of $\mathrm{\mathbf{U}}_{r}$, $\mathrm{\mathbf{V}}_{r}$ and $\mathrm{\mathbf{W}}_{r}$ belong to $q\mathbb{Z}$ for some $q\in\mathbb{R}$.

Discrete decompositions are favorable for two reasons. The first reason concerns the exactness of the decomposition: if $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ is computed with numerical methods, then this will give a decomposition of $\Phi_{m,p,n}$ only up to some finite accuracy (and also up to machine precision, due to floating-point arithmetic computations). Hence, the matrix multiplication algorithm obtained from this decomposition will compute the product of $\mathrm{\mathbf{A}}\in\mathbb{R}^{m\times p}$ and $\mathrm{\mathbf{B}}\in\mathbb{R}^{p\times n}$ with a small error, even in exact arithmetic. This is not advisable because this error will accumulate when we will apply the algorithm in a recursive way to compute the product of general $N\times N$ matrices [2, Chapter 15]. These limited-accuracy issues can be overcome if we know a priori that $\mathrm{\mathbf{U}}_{r}$, $\mathrm{\mathbf{V}}_{r}$ and $\mathrm{\mathbf{W}}_{r}$ have their entries in a known discrete set.

The second reason to favor discrete decompositions is that the obtained algorithm for matrix multiplication will have better stability and computational cost properties. Indeed, if the entries of $\mathrm{\mathbf{U}}_{r}$, $\mathrm{\mathbf{V}}_{r}$ and $\mathrm{\mathbf{W}}_{r}$ belong to $q\mathbb{Z}$, then it is not hard to show (we do not go into the details) that, modulo some pre- and post multiplication of $\mathrm{\mathbf{A}}\in\mathbb{R}^{m\times p}$ and $\mathrm{\mathbf{B}}\in\mathbb{R}^{p\times n}$ by $q$, the product $\mathrm{\mathbf{A}}\mathrm{\mathbf{B}}$ can be computed using only additions and multiplications of the entries of $\mathrm{\mathbf{A}}$ and $\mathrm{\mathbf{B}}$ by integers. (For example, in Strassen’s algorithm, $f_{r}(\mathrm{\mathbf{A}})$ [resp. $g_{r}(B)$] can be obtained using only additions and subtractions of the entries of $\mathrm{\mathbf{A}}$ [resp. $\mathrm{\mathbf{B}}$].) Multiplication by integers is more rapid and stable than multiplication by arbitrary floating-point numbers (for instance, multiplication by a power of $2$ is equivalent to changing the exponent in the floating point representation). For more detailed information on the forward normwise error induced by a fast matrix multiplication algorithm, we refer the interested reader to [13].

3 Invariance transformations

The main goal of this paper is to study relations between decompositions of a given matrix multiplication tensor. We describe three types of operations that transform an $F$-PD of a matrix multiplication tensor into another $F$-PD of the same tensor. These transformations will be referred to as invariance transformations.

Proposition 3.1** (Invariance transformations).**

Let $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ be an $F$-PD of $\Phi_{m,p,n}$. The following transformations produce matrices $\mathrm{\mathbf{U}}^{\prime}_{r}$, $\mathrm{\mathbf{V}}^{\prime}_{r}$ and $\mathrm{\mathbf{W}}^{\prime}_{r}$ ($1\leq r\leq F$) such that $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ is also an $F$-PD of $\Phi_{m,p,n}$:

• •

Permutation transformations: let $\sigma\in\mathfrak{S}_{F}$ (where $\mathfrak{S}_{F}$ is the set of permutations of $\{1,\ldots,F\}$) and define

[TABLE]

• •

Scaling transformations: choose coefficients $\lambda_{r},\mu_{r},\nu_{r}\in\mathbb{R}$ such that $\lambda_{r}\mu_{r}\nu_{r}=1$ for each $1\leq r\leq F$ and define

[TABLE]

• •

Trace transformations: let $\mathrm{\mathbf{P}}\in\mathrm{GL}(m)$, $\mathrm{\mathbf{Q}}\in\mathrm{GL}(p)$ and $\mathrm{\mathbf{R}}\in\mathrm{GL}(n)$ (where $\mathrm{GL}(h)$ denotes the set of invertible $h\times h$ matrices), and define

[TABLE]

The first two classes of invariance transformations (permutations and scaling) provide invariance transformations for the decompositions of any tensors. However, the third class (trace transformations) is somehow specific to matrix multiplication tensors, as it originates from the invariance of the trace operator (see proof below); hence the name “trace transformations”.

Proof of Proposition 3.1.

(See, e.g., [22]). The invariance of the permutation and scaling transformations is straightforward. For the trace transformations, let $f^{\prime}_{r}\in(\mathbb{R}^{m\times p})^{*}$, $g^{\prime}_{r}\in(\mathbb{R}^{p\times n})^{*}$ and $\Phi^{\prime}\in\mathrm{Bil}(\mathbb{R}^{m\times p},\mathbb{R}^{p\times n};\mathbb{R}^{m\times n})$ be given by (2) and (1) with $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$. Then

[TABLE]

where we have used the invariance property of the trace operator with respect to cyclic permutations. Similarly, $g^{\prime}_{r}(\mathrm{\mathbf{B}})=g_{r}(\mathrm{\mathbf{Q}}\mathrm{\mathbf{B}}\mathrm{\mathbf{R}}^{-1})$. It follows that

[TABLE]

Thus $\Phi^{\prime}=\Phi_{m,p,n}$ showing that $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ is an $F$-PD of $\Phi_{m,p,n}$. ∎

Invariance transformations define an equivalence relation on the set of $F$-PDs of a given matrix multiplication tensor. For fixed $m,p,n$ and $F$, two polyadic decompositions $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ of $\Phi_{m,p,n}$ are equivalent if there exist permutation, scaling and/or trace transformations that allow one to transform $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ into $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$. We will also say that $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ are permutation-equivalent if there exists a permutation transformation allowing us to transform $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ into $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$. Similarly, we define the notion of (scaling+trace)-equivalence.

We have just seen that invariance transformations can be used to produce many different decompositions (i.e., many fast matrix multiplication algorithms) from a given one. This raises the following questions that we will address in this paper:

How many inequivalent polyadic decompositions does $\Phi_{m,p,n}$ admit? In other words, how many essentially different fast matrix multiplication algorithms are there for the multiplication of $m\times p$ matrices by $p\times n$ matrices? We will tackle this question in Sections 4 and 6.3.

Starting from a given algorithm for the multiplication of $m\times p$ matrices by $p\times n$ matrices, can we obtain with invariance transformations another algorithm with better performance (e.g., in terms of stability and efficiency)? We will tackle this question in Sections 5 and 6.2.

Remark 3.1*.*

At first sight, it might look like the scaling transformations act as a particular case of the trace transformations with

[TABLE]

for example. In fact, this is not the case since the above $\mathrm{\mathbf{P}},\mathrm{\mathbf{Q}},\mathrm{\mathbf{R}}$ will rescale all the matrices $\mathrm{\mathbf{U}}_{r},\mathrm{\mathbf{V}}_{r},\mathrm{\mathbf{W}}_{r}$ with the same coefficients $\lambda,\mu,\nu$ (provided $\lambda\mu\nu=1$) while the scaling transformations admit different coefficients $\lambda_{1},\ldots,\lambda_{F}$, $\mu_{1},\ldots,\mu_{F}$ and $\nu_{1},\ldots,\nu_{F}$.  $\triangleleft$

4 An algorithm for checking equivalence

In this section, we present an algorithm for deciding whether two $F$-PDs of a matrix multiplication tensor are equivalent. Under mild assumptions on the input $F$-PDs, the algorithm will either return the permutation, scaling and trace transformations that allow one to connect both $F$-PDs or conclude that the two $F$-PDs are not equivalent to each other. The working assumptions were satisfied for $100\%$ of the samples on which we performed numerical experiments (see Section 6.3), motivating the qualifier “mild” assumptions.

We start this section by introducing the concept of clustering number of a matrix. This number can be computed efficiently and is used in the assumption to guarantee proper working of the algorithm.

4.1 The clustering number of a matrix

Let $\mathrm{\mathbf{A}}$ be an $m\times n$ matrix. Let $\{U_{1},\ldots,U_{S}\}$ be a family of linearly independent subspaces of $\mathbb{R}^{m}$ (i.e., $u_{1}+\cdots+u_{S}=0$ with $u_{i}\in U_{i}$ for each $1\leq i\leq S$ implies $u_{i}=0$ for every $1\leq i\leq S$). If each column of $\mathrm{\mathbf{A}}$ belongs to some $U_{i}$ (in fact, except the case where the column contains only zeros, it may belong to at most one $U_{i}$ since they are linearly independent), then we say that $\{U_{1},\ldots,U_{S}\}$ is a cover of $\mathrm{\mathbf{A}}$. The largest integer $S^{*}$ such that there exists a cover of $\mathrm{\mathbf{A}}$ with $S^{*}$ linearly independent subspaces is called the clustering number of $\mathrm{\mathbf{A}}$ and is denoted by $\mathrm{cl}_{\oplus}(A)=S^{*}$ (the choice of the symbol $\oplus$ comes from the fact that if $\{U_{1},\ldots,U_{S}\}$ is a cover of $A$ with $S=S^{*}$ then $U_{1}\oplus\cdots\oplus U_{S}=\mathbb{R}^{m}$).

Example 4.1*.*

Let $U_{1}=\mathrm{colspan}(\mathrm{\mathbf{A}})$ and suppose that the rank of $\mathrm{\mathbf{A}}$ is $r$. Then it is easy to build a cover $\{U_{1},U^{\prime}_{1},\ldots,U^{\prime}_{m-r}\}$ where the $U^{\prime}_{i}$’s are one-dimensional subspaces. Hence, we conclude that $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})\geq m+1-\mathrm{rank}(\mathrm{\mathbf{A}})$.  $\triangleleft$

Suppose that $\mathrm{\mathbf{A}}$ has full row-rank, i.e., $\mathrm{rank}(\mathrm{\mathbf{A}})=m$. We give a characterization of the clustering number of $\mathrm{\mathbf{A}}$ in terms of the connected components of a graph. Denote by $\mathrm{\mathbf{a}}_{1},\ldots,\mathrm{\mathbf{a}}_{n}$ the columns of $\mathrm{\mathbf{A}}$. Without loss of generality, we may assume that the first $m$ columns of $\mathrm{\mathbf{A}}$ span $\mathbb{R}^{m}$. Let $\mathrm{\mathbf{A}}^{\prime}=[\mathrm{\mathbf{a}}_{1},\ldots,\mathrm{\mathbf{a}}_{m}]$ and $\mathrm{\mathbf{A}}^{\prime\prime}=[\mathrm{\mathbf{a}}_{m+1},\ldots,\mathrm{\mathbf{a}}_{n}]$.

Define the undirected graph G as follows. The integers $\{1,\ldots,m\}$ are the nodes of G and for each column $\mathrm{\mathbf{a}}_{j}$ of $\mathrm{\mathbf{A}}^{\prime\prime}$, let $\mathrm{\mathbf{q}}_{j}=(q_{j,1},\ldots,q_{j,m})^{\top}$ be its coordinates in the basis defined by $\mathrm{\mathbf{A}}^{\prime}$, i.e., $\mathrm{\mathbf{q}}_{j}=(\mathrm{\mathbf{A}}^{\prime})^{-1}\mathrm{\mathbf{a}}_{j}$. For each $i_{1},i_{2}\in\{1,\ldots,m\}$, draw an edge between the nodes $i_{1}$ and $i_{2}$ if and only if there is a column $\mathrm{\mathbf{a}}_{j}$ of $\mathrm{\mathbf{A}}^{\prime\prime}$ such that its coordinate vector has a nonzero component in both $\mathrm{\mathbf{a}}_{i_{1}}$ and $\mathrm{\mathbf{a}}_{i_{2}}$, i.e., if $q_{j,i_{1}}\neq 0$ and $q_{j,i_{2}}\neq 0$ (clearly, there might be multiple edges and also loops). Moreover, each edge receives a label: this label is simply $j$, the index of the column of $\mathrm{\mathbf{A}}^{\prime\prime}$ that led to this edge. An example is represented in Figure 1. This construction allows us to state the following lemma:

Proposition 4.1**.**

Let $\mathrm{\mathbf{A}}$ and G be defined as above. Then the clustering number of $\mathrm{\mathbf{A}}$ is equal to the number of connected components of G.

Proof.

Let $\{\mbox{G}_{1},\ldots,\mbox{G}_{T}\}$ be the connected components of G. First, we show that $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})\geq T$. For each $1\leq t\leq T$, let $I_{t}$ be the set of all column indices involved in the nodes of $\mbox{G}_{t}$. (For example, considering the graph in Figure 1, we would have $I_{1}=\{1,2\}$ and $I_{2}=\{3\}$.) For each $t$, define the subspace $U_{t}$ as the subspace spanned by the columns $\mathrm{\mathbf{a}}_{i}$ with $i\in I_{t}$. By hypothesis on $\mathrm{\mathbf{A}}^{\prime}$ having full rank, the subspaces $U_{t}$ satisfy $\mathbb{R}^{m}=U_{1}\oplus\cdots\oplus U_{T}$. We have to show that each column $\mathrm{\mathbf{a}}_{j}$ of $\mathrm{\mathbf{A}}^{\prime\prime}$ belongs to some $U_{t}$.

Therefore, we show that the label $j$ appears in the edges of at most one component $\mbox{G}_{t}$. Indeed, if $j$ appears in $\mbox{G}_{t_{1}}$ and $\mbox{G}_{t_{2}}$, then $\mathrm{\mathbf{a}}_{j}$ has a nonzero component in at least one node $i_{1}$ of $\mbox{G}_{t_{1}}$ and one node $i_{2}$ of $\mbox{G}_{t_{2}}$. Hence, there must be an edge between $i_{1}$ and $i_{2}$ and thus $\mbox{G}_{t_{1}}$ and $\mbox{G}_{t_{2}}$ are connected, a contradiction. Thus $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})\geq T$.

To show that $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})\leq T$, let $S=\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})$ and let $\{U_{1},\ldots,U_{S}\}$ be a cover of $\mathbb{R}^{m}$. For each $1\leq s\leq S$, let $I_{s}$ be the set of the indices of the columns of $\mathrm{\mathbf{A}}^{\prime}$ belonging to $U_{s}$. Since $\mathrm{\mathbf{A}}^{\prime}$ has full rank, it is clear that $U_{s}=\mathrm{span}(\{\mathrm{\mathbf{a}}_{i}\}_{i\in I_{s}})$. We show that $\{I_{1},\ldots,I_{S}\}$ defines connected components of G. Indeed, if there is an edge, say with label $j$, between nodes $i_{1}\in I_{s_{1}}$ and $i_{2}\in I_{s_{2}}$, then $\mathrm{\mathbf{a}}_{j}$ has nonzero components in $\mathrm{\mathbf{a}}_{i_{1}}$ and in $\mathrm{\mathbf{a}}_{i_{2}}$ and thus it belongs to the subspaces $U_{s_{1}}$ and $U_{s_{2}}$. However, by definition, the subspaces $U_{s}$ are linearly independent. Hence, we must have $s_{1}=s_{2}$. We conclude that there are at least $S$ connected components in G and thus $T\geq\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})$. ∎

Proposition 4.1 above gives an efficient way to compute the clustering number of a matrix. We describe below another way to efficiently compute the clustering number of a matrix using linear algebra only (and that will be useful in the proof of Lemma 4.5 below). To do this, let $\mathrm{\mathbf{A}}$ be an $m\times n$ matrix and consider the following linear system:

[TABLE]

with variables $\mathrm{\mathbf{M}}\in\mathbb{R}^{m\times m}$ and $\xi_{[n]}=(\xi_{1},\ldots,\xi_{n})\in\mathbb{R}^{n}$. Clearly the system (3) is homogeneous. Let us denote by $\mathscr{S}$ the vector space of solutions $(\mathrm{\mathbf{M}},\xi_{[n]})$ to (3).

Lemma 4.2**.**

Let $\mathrm{\mathbf{A}}$ have full row-rank and no zero columns and let $\mathscr{S}$ be defined as above. Then $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})=\mathrm{dim}(\mathscr{S})$.

Proof.

First note that if $\mathrm{\mathbf{A}}\in\mathbb{R}^{m\times n}$ has full row-rank then $m\leq n$. Let $\mathrm{\mathbf{A}}^{\prime}$ and $\mathrm{\mathbf{A}}^{\prime\prime}$ be defined as previously. Without loss of generality, we may again assume that $\mathrm{\mathbf{A}}^{\prime}$ has full rank. Let $\xi_{[n]}=(\xi_{1},\ldots,\xi_{n})$ be fixed. Then we have

[TABLE]

whence $\mathrm{\mathbf{M}}$ is completely determined by $\xi_{[n]}$. Reversely, if $\mathrm{\mathbf{M}}$ is known, then every $\xi_{1},\ldots,\xi_{n}$ are also determined since $\mathrm{\mathbf{A}}$ has no zero columns. Hence, we only have to compute the number of degrees of freedom in $\xi_{[n]}$ to compute the dimension of $\mathscr{S}$.

Let G be the graph associated to $\mathrm{\mathbf{A}}$. Suppose that $i_{1}$ and $i_{2}$ are two nodes that are adjacent to each other with an edge labeled by $j$. Then, we have

[TABLE]

Since $\mathrm{\mathbf{a}}_{1},\ldots,\mathrm{\mathbf{a}}_{m}$ are linearly independent and $q_{j,i_{1}}\neq 0$ and $q_{j,i_{2}}\neq 0$, the only solution provides that $\xi_{i_{1}}=\xi_{i_{2}}=\xi_{j}$. We conclude that if two nodes $i_{1},i_{2}$ belong to the same connected component $\mbox{G}_{t}$, then $\xi_{i_{1}}=\xi_{i_{2}}$. Hence, using Lemma 4.1, the dimension of $\mathscr{S}$ is lower than or equal to $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})$.

On the other hand, let $S=\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})$ and let $\{U_{1},\ldots,U_{S}\}$ be a cover of $\mathrm{\mathbf{A}}$. Then for each $1\leq s\leq S$, let $\mathrm{\mathbf{M}}_{s}$ be the projection on $U_{s}$, i.e.,

[TABLE]

Let $(\eta_{1},\ldots,\eta_{S})$ be a fixed vector and define $\mathrm{\mathbf{M}}=\sum_{s=1}^{S}\eta_{s}\mathrm{\mathbf{M}}_{s}$. For each $1\leq i\leq n$, let $\xi_{i}=\eta_{s_{i}}$ where $s_{i}$ is the unique index $1\leq s_{i}\leq S$ such that $\mathrm{\mathbf{a}}_{i}\in U_{s_{i}}$. Then $(\mathrm{\mathbf{M}},\xi_{[n]})$ is a solution of (3). Hence, the dimension of $\mathscr{S}$ is at least $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})$. ∎

We are now able to state and prove the main theorem for this subsection:

Theorem 4.3**.**

Let $\mathrm{\mathbf{A}}$ be an $m\times n$ matrix with rank $r$ and let $Z$ be the number of zero columns in $\mathrm{\mathbf{A}}$. Then the clustering number of $\mathrm{\mathbf{A}}$ and the dimension of the solution space $\mathscr{S}$ of (3) satisfy

[TABLE]

Proof.

First we suppose that there are no zero columns in $\mathrm{\mathbf{A}}$. Let $\mathrm{\mathbf{X}}\in\mathrm{GL}(m)$ and observe that $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{X}}\mathrm{\mathbf{A}})=\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})$. Consider the linear system

[TABLE]

and note that $(\mathrm{\mathbf{M}},\xi_{[n]})$ is a solution of (4) if and only if $(\mathrm{\mathbf{X}}^{-1}\mathrm{\mathbf{M}}\mathrm{\mathbf{X}},\xi_{[n]})$ is a solution of (3). Hence, taking an appropriate matrix $\mathrm{\mathbf{X}}$, we may assume without loss of generality that the last $m-r$ rows of $\mathrm{\mathbf{A}}$ are zero.

Let $\tilde{\mathrm{\mathbf{A}}}$ be the $r\times n$ matrix consisting of the first $r$ rows of $\mathrm{\mathbf{A}}$. Then $\tilde{\mathrm{\mathbf{A}}}$ has full row-rank and it is easy to check that $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})=\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{A}}})+m-r$. The matrix $\mathrm{\mathbf{M}}$ may be partitioned into the following blocks:

[TABLE]

It is clear that $(\mathrm{\mathbf{M}},\xi_{[n]})$ is a solution of (3) if and only if

[TABLE]

From Lemma 4.2 and the fact that $\tilde{\mathrm{\mathbf{A}}}$ has full row-rank, the solutions $(\mathrm{\mathbf{M}}_{1},\xi_{[n]},\mathrm{\mathbf{M}}_{2})$ of (5) form a vector space with dimension $\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{A}}})$. On the other hand, there are no constraints on $\mathrm{\mathbf{M}}_{3}\in\mathbb{R}^{m\times(m-r)}$. This proves the assertion when $Z=0$.

Finally, observe that appending a zero column to $\mathrm{\mathbf{A}}$ does not change its clustering number and also does not change the space of solutions $\mathrm{\mathbf{M}}$ of (3). The only thing that changes is that the coefficient $\xi_{n+1}$ affected to this zero column might take any value. Hence, the dimension of $\mathscr{S}$ is increased by one. This concludes the proof of the theorem. ∎

Remark 4.1*.*

For the interested reader, let us mention that the clustering number has an interpretation in terms of matroids: considering $\mathrm{\mathbf{A}}$ as a linear matroid with ground set given by the columns of $\mathrm{\mathbf{A}}$ [29, Chapter 39], then we can show that the clustering number of $\mathrm{\mathbf{A}}$ is in fact equal to the number of connected components (in the matroid sense [29, Chapter 39]) of the matroid $\mathrm{\mathbf{A}}$. Dedicated softwares exist to compute the connected components of a matroid. See, e.g., [30]. In fact, in the case of a linear matroid, the implementation in [30] is equivalent to dynamically computing the connected components of the graph G in Proposition 4.1.  $\triangleleft$

4.2 Computation of the scaling and trace transformations

Let $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ be two $F$-PDs of a matrix multiplication tensor $\Phi_{m,p,n}$. We would like to know whether they are equivalent (see Section 3) and, if they are, to compute the invariance transformations connecting $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ to $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$. At first, we assume that the permutation transformation is given and we focus on the computation of the scaling and trace transformations. (We will see in the next subsection how we can compute this permutation transformation without trying all permutations of $\{1,\ldots,F\}$.) Under mild assumptions on $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$, we will see how to do this computation using linear algebra only.

First, we make two important comments. In the sequel, we will always assume that the rank-$1$ terms555$f_{r}\otimes g_{r}$ denotes the tensor product of functions $f_{r}$ and $g_{r}$, i.e., $(f_{r}\otimes g_{r})(u,v)=f_{r}(u)g_{r}(v)$. $(f_{r}\otimes g_{r})w_{r}$, $1\leq r\leq F$, in the $F$-PDs (1) are linearly independent. Indeed, if one of the rank-$1$ terms $(f_{r}\otimes g_{r})w_{r}$ can be decomposed as a linear combination of the other rank-$1$ terms, then it is easy to build a polyadic decomposition (1) of $\Phi_{m,p,n}$ with $F-1$ terms (which would be extremely lucky and never happened in the numerical experiments we performed; see Section 6).666Even in this eventuality, if the decompositions $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ are equivalent and if the rank-$1$ terms of $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ are linearly dependent, then the rank-$1$ terms of $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ are also linearly dependent. Moreover, the $F_{*}$-term polyadic decompositions $(\mathrm{\mathbf{U}}_{[F^{*}]},\mathrm{\mathbf{V}}_{[F^{*}]},\mathrm{\mathbf{W}}_{[F^{*}]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F^{*}]},{\mathrm{\mathbf{V}}^{\prime}}_{[F^{*}]},{\mathrm{\mathbf{W}}^{\prime}}_{[F^{*}]})$, $1\leq F_{*}<F$, obtained by removing in $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ the linearly dependent rank-$1$ terms and the corresponding terms in $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$, are equivalent and contain both linearly independent rank-$1$ terms. Hence, modulo a little extra work (due to the non-uniqueness in the choice of linearly dependent terms to remove), we may always reduce to the case where $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ contains only linearly independent rank-$1$ terms.

The second comment is summarized in the following theorem:

Theorem 4.4**.**

Let $\Phi_{m,p,n}$ be a matrix multiplication tensor and let $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ be an $F$-PD of $\Phi_{m,p,n}$. Let $I\subseteq\{1,\ldots,F\}$ be a subset of indices with $\lvert I\rvert+n\geq F+1$. Then the family $\{\mathrm{\mathbf{U}}_{r}\}_{r\in I}$ fully spans $\mathbb{R}^{p\times m}$.

Proof.

Suppose, on the contrary, that there exists a subset $I\subseteq\{1,\ldots,F\}$ with size $\ell$ such that $\ell+n\geq F+1$ and $\mathrm{span}(\{\mathrm{\mathbf{U}}_{r}\}_{r\in I})\neq\mathbb{R}^{p\times m}$. Then there exists $\mathrm{\mathbf{A}}_{*}\in\mathbb{R}^{m\times p}\setminus\{0\}$ such that $\mathrm{trace}(\mathrm{\mathbf{U}}_{r}\mathrm{\mathbf{A}}_{*})=0$ for each $r\in I$. Denote by $Y$ the vector space of $p\times n$ matrices $\mathrm{\mathbf{B}}$ such that $\mathrm{trace}(\mathrm{\mathbf{V}}_{r}\mathrm{\mathbf{B}})=0$ for each $r\in F\setminus I$. Since $\lvert F\setminus I\rvert=F-\ell\leq n-1$, we have that $\mathrm{dim}(Y)\geq pn-n+1$.

Now define $Z$ as the vector space of all matrices $\mathrm{\mathbf{Z}}=(\mathrm{\mathbf{C}}\mathrm{\mathbf{A}}_{*})^{\top}\in\mathbb{R}^{p\times n}$ for some $\mathrm{\mathbf{C}}\in\mathbb{R}^{n\times m}$. Since $\mathrm{\mathbf{A}}_{*}\neq 0$, the dimension of $Z$ is at least $n$. Now let $\mathrm{\mathbf{B}}\in Y$ and $\mathrm{\mathbf{Z}}=(\mathrm{\mathbf{C}}\mathrm{\mathbf{A}}_{*})^{\top}\in Z$ and observe that

[TABLE]

From the definitions of $\mathrm{\mathbf{A}}_{*}$ and $Y$, we conclude that $\mathrm{trace}(\mathrm{\mathbf{B}}\mathrm{\mathbf{Z}}^{\top})=0$. Thus $Y\cap Z=\{0\}$, so $\mathrm{dim}(Y)+\mathrm{dim}(Z)=pn-n+1+n>pn$. This contradicts $Y,Z\subseteq\mathbb{R}^{p\times n}$. ∎

Remark 4.2*.*

Note that Theorem 4.4 applies, mutatis mutandis, to $\mathrm{\mathbf{V}}_{[F]}$ and $\mathrm{\mathbf{W}}_{[F]}$. To see this, it suffices to observe that if $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ is an $F$-PD of $\Phi_{m,p,n}$ then $(\mathrm{\mathbf{W}}_{[F]},\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]})$ and $(\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]},\mathrm{\mathbf{U}}_{[F]})$ provide $F$-PDs of $\Phi_{n,m,p}$ and $\Phi_{p,n,m}$ respectively.  $\triangleleft$

Among other conclusions of this theorem, we get that $\mathrm{span}(\{\mathrm{\mathbf{U}}_{r}\}_{1\leq r\leq F})=\mathbb{R}^{p\times m}$. Indeed, it suffices to apply Theorem 4.4 with $I=\{1,\ldots,F\}$. Similar conclusions hold for $\mathrm{\mathbf{V}}_{[F]}$ and $\mathrm{\mathbf{W}}_{[F]}$.

We now present the algorithm to compute the scaling and trace transformations between $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ or conclude that no such transformations exist. To simplify the notation, it will be useful to consider $\mathrm{\mathbf{U}}_{r}$, $\mathrm{\mathbf{V}}_{r}$ and $\mathrm{\mathbf{W}}_{r}$ as column vectors and gather them into matrices. Therefore, we define

[TABLE]

where $\mathrm{vec}(\cdot)$ is the vectorization (column stacking) operator. Similarly, we define $\tilde{\mathrm{\mathbf{V}}}\in\mathbb{R}^{np\times F}$ and $\tilde{\mathrm{\mathbf{W}}}\in\mathbb{R}^{mn\times F}$, and also $\tilde{\mathrm{\mathbf{U}}}^{\prime}\in\mathbb{R}^{pm\times F}$, $\tilde{\mathrm{\mathbf{V}}}^{\prime}\in\mathbb{R}^{np\times F}$ and $\tilde{\mathrm{\mathbf{W}}}^{\prime}\in\mathbb{R}^{mn\times F}$. The algorithm is guaranteed to work if we make the following assumption on $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$:

Assumption 4.1**.**

Let $\tilde{\mathrm{\mathbf{U}}}$, $\tilde{\mathrm{\mathbf{V}}}$ and $\tilde{\mathrm{\mathbf{W}}}$ be defined as above. We assume that either $\tilde{\mathrm{\mathbf{U}}}$, $\tilde{\mathrm{\mathbf{V}}}$ or $\tilde{\mathrm{\mathbf{W}}}$ has clustering number equal to one.

Remark 4.3*.*

It is not difficult to see that if $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ are equivalent, then $\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{U}}})=\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{U}}}^{\prime})$, $\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{V}}})=\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{V}}}^{\prime})$ and $\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{W}}})=\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{W}}}^{\prime})$. Thus the clustering numbers (which can be efficiently computed) already offer us a way to eliminate $F$-PDs that are not equivalent. Therefrom, Assumption 4.1 can be rephrased (without loss of generality) as follows: either $\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{U}}})=\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{U}}}^{\prime})=1$, or $\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{V}}})=\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{V}}}^{\prime})=1$, or $\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{W}}})=\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{W}}}^{\prime})=1$.  $\triangleleft$

We will see in Section 6.3 that Assumption 4.1 is satisfied for $100\%$ of the randomly computed samples on which we have performed numerical experiments. The goal of the algorithm presented in this subsection is to compute matrices $\mathrm{\mathbf{P}}\in\mathrm{GL}(m)$, $\mathrm{\mathbf{Q}}\in\mathrm{GL}(p)$ and $\mathrm{\mathbf{R}}\in\mathrm{GL}(n)$, and scaling coefficients $\lambda_{1},\ldots,\lambda_{F}$, $\mu_{1},\ldots,\mu_{F}$ and $\nu_{1},\ldots,\nu_{F}$ such that $\lambda_{r}\mu_{r}\nu_{r}=1$ and

[TABLE]

for every $1\leq r\leq F$. The above conditions are nonlinear in $\mathrm{\mathbf{P}}$, $\mathrm{\mathbf{Q}}$, $\mathrm{\mathbf{R}}$, $\lambda_{r}$, $\mu_{r}$ and $\nu_{r}$. However, relying on the assumptions on $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$, these conditions can be reduced to linear matrix equations.

First of all, we show that the requirement $\lambda_{r}\mu_{r}\nu_{r}=1$ can be dropped. Indeed, suppose that $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ satisfy (7), and let $f^{\prime}_{r}$, $g^{\prime}_{r}$ and $w^{\prime}_{r}$ be given by (2) with $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$. Also let $\Phi^{\prime\prime}$ be given by (1) and (2) with $\mathrm{\mathbf{U}}^{\prime\prime}_{r}=\mathrm{\mathbf{Q}}^{-1}\mathrm{\mathbf{U}}_{r}\mathrm{\mathbf{P}}$, $\mathrm{\mathbf{V}}^{\prime\prime}_{r}=\mathrm{\mathbf{R}}^{-1}\mathrm{\mathbf{V}}_{r}\mathrm{\mathbf{Q}}$, $\mathrm{\mathbf{W}}^{\prime\prime}_{r}=\mathrm{\mathbf{P}}^{-1}\mathrm{\mathbf{W}}_{r}\mathrm{\mathbf{R}}$. Then $\Phi^{\prime\prime}=\Phi_{m,p,n}$ (trace transformations) and

[TABLE]

From the linear independence assumption on the rank-$1$ terms $(f^{\prime}_{r}\otimes g^{\prime}_{r})w^{\prime}_{r}$, $1\leq r\leq F$, we conclude that $\lambda_{r}\mu_{r}\nu_{r}=1$ is trivially satisfied if (7) holds.

According to Assumption 4.1, we assume for the rest of this subsection that $\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{U}}})=\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{U}}}^{\prime})=1$. We denote by $\mathrm{\mathbf{A}}\otimes\mathrm{\mathbf{B}}$ the Kronecker product of two matrices $\mathrm{\mathbf{A}}$ and $\mathrm{\mathbf{B}}$, and we will use the following property of the vectorization operator:

[TABLE]

Then the first equation of (7) is equivalent to

[TABLE]

Considering $\mathrm{\mathbf{P}}^{\top}\!\otimes\mathrm{\mathbf{Q}}^{-1}$ as a single matrix $\mathrm{\mathbf{M}}\in\mathbb{R}^{pm\times pm}$, (8) becomes

[TABLE]

which is linear in $\mathrm{\mathbf{M}}$ and $\lambda_{[F]}=(\lambda_{1},\ldots,\lambda_{F})$. The fact that no unwanted solutions are created by this linearization is shown in the following developments.

Let $\mathrm{\mathbf{A}}$ and $\mathrm{\mathbf{A}}^{\prime}$ be two $m\times n$ matrices with full row-rank, containing no zero columns and with $\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}})=\mathrm{cl}_{\oplus}(\mathrm{\mathbf{A}}^{\prime})=1$. Then consider the linear system

[TABLE]

with variables $\mathrm{\mathbf{M}}\in\mathbb{R}^{m\times m}$ and $\xi_{[n]}=(\xi_{1},\ldots,\xi_{n})\in\mathbb{R}^{n}$. This problem is close to problem (3) except that we allow $\mathrm{\mathbf{A}}\neq\mathrm{\mathbf{A}}^{\prime}$. Let $\mathscr{S}$ be the vector space of $(\mathrm{\mathbf{M}},\xi_{[n]})$ that are solutions of (10).

Lemma 4.5**.**

Let $\mathscr{S}$ be defined as above. If $\mathscr{S}$ contains a solution $(\mathrm{\mathbf{M}},\xi_{[n]})$ such that $\xi_{i}\neq 0$ for every $1\leq i\leq n$, then $\mathrm{dim}(\mathscr{S})=1$.

Proof.

Let $(\mathrm{\mathbf{M}},\xi_{[n]})$ be a solution of (10) with $\xi_{i}\neq 0$ for every $1\leq i\leq n$. We have assumed that $\mathrm{\mathbf{A}}^{\prime}$ has full row-rank and thus $\mathrm{\mathbf{A}}^{\prime}\,\mathrm{diag}(\xi_{1},\ldots,\xi_{n})$ has full row-rank as well. Hence, $\mathrm{\mathbf{M}}$ must be invertible. In a similar way as in the proof of Lemma 4.2, we may assume without loss of generality that the first $m$ columns of $\mathrm{\mathbf{A}}$ span $\mathbb{R}^{m}$. Hence, the first $m$ columns of $\mathrm{\mathbf{A}}^{\prime}$ span $\mathbb{R}^{m}$ too. We conclude the proof with a similar reasoning as for the first part of the proof of Lemma 4.2. ∎

Hence, two cases can happen when solving (8): (i) either the linearized system (9) admits no solutions with $\lambda_{r}\neq 0$ for every $1\leq r\leq F$; in this case, we conclude that the two $F$-PDs are not (scaling+trace)-equivalent; or (ii) the solution space $\mathscr{S}$ of (9) is one-dimensional and thus taking an arbitrary nonzero $(\mathrm{\mathbf{M}},\lambda_{[F]})\in\mathscr{S}$, it is easy to check whether $\mathrm{\mathbf{M}}$ has the form $\mathrm{\mathbf{M}}=\mathrm{\mathbf{P}}^{\top}\!\otimes\mathrm{\mathbf{Q}}^{-1}$ for some $\mathrm{\mathbf{P}}\in\mathrm{GL}(m)$ and $\mathrm{\mathbf{Q}}\in\mathrm{GL}(p)$. If the latter does not hold, then the two $F$-PDs are not (scaling+trace)-equivalent. Otherwise, $\mathrm{\mathbf{P}}$ and $\mathrm{\mathbf{Q}}$ are the unique (up to a scalar multiplication) matrices involved in the invariance transformations (7).

Now that we have determined $\mathrm{\mathbf{P}}$ and $\mathrm{\mathbf{Q}}$, we consider the following linear system:

[TABLE]

where the unknowns are $\mathrm{\mathbf{R}}\in\mathbb{R}^{n\times n}$ and $\tilde{\mu}_{[F]},\tilde{\nu}_{[F]}\in\mathbb{R}^{F}$ for $1\leq r\leq F$. If (11) admits no solutions $(\mathrm{\mathbf{R}},\tilde{\mu}_{[F]},\tilde{\nu}_{[F]})$ with $\tilde{\mu}_{r}\neq 0$ and $\tilde{\nu}_{r}\neq 0$ for every $1\leq r\leq F$, then we conclude that $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ are not (scaling+trace)-equivalent. On the other hand, if $\tilde{\mu}_{r}\neq 0$ and $\tilde{\nu}_{r}\neq 0$ for every $1\leq r\leq F$, then $\mathrm{\mathbf{R}}$ is invertible because $\mathrm{span}(\{\mathrm{\mathbf{V}}_{r}\mathrm{\mathbf{Q}}\}_{1\leq r\leq F})=\mathbb{R}^{n\times p}$ (see Remark 4.2). The $6$-tuple $(\mathrm{\mathbf{P}},\mathrm{\mathbf{Q}},\mathrm{\mathbf{R}},\lambda_{[F]},\mu_{[F]},\nu_{[F]})$ with $\mu_{r}=\tilde{\mu}_{r}^{-1}$ and $\nu_{r}=\tilde{\nu}_{r}^{-1}$ then provides a solution to the (scaling+trace)-equivalence problem (7).

4.3 Computation of the permutation transformation

In the previous subsection, we have described a procedure to compute the scaling and trace transformations connecting two $F$-PDs $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ or conclude that no such transformations exist. The equivalence of $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ can then be decided in finite time by trying every permutation $\sigma\in\mathfrak{S}_{F}$ and testing the (scaling+trace)-equivalence of $\sigma((\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]}))$777where $\sigma((\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]}))=(\sigma(\mathrm{\mathbf{U}}_{[F]}),\sigma(\mathrm{\mathbf{V}}_{[F]}),\sigma(\mathrm{\mathbf{W}}_{[F]}))$ and $\sigma(\mathrm{\mathbf{X}}_{[F]})$ is the permuted $F$-uple $(\mathrm{\mathbf{X}}_{\sigma(1)},\ldots,\mathrm{\mathbf{X}}_{\sigma(F)})$. and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$. Due to the combinatorial growth of $\lvert\mathfrak{S}_{F}\rvert$, an exhaustive exploration of $\mathfrak{S}_{F}$ is generally not feasible in practice. In this section, we explain how to efficiently decide whether the two $F$-PDs are equivalent without trying all permutations $\sigma\in\mathfrak{S}_{F}$.

Definition 4.1**.**

Let $\mathrm{\mathbf{A}}_{[m]}=(\mathrm{\mathbf{A}}_{1},\ldots,\mathrm{\mathbf{A}}_{m})$ and $\mathrm{\mathbf{B}}_{[m]}=(\mathrm{\mathbf{B}}_{1},\ldots,\mathrm{\mathbf{B}}_{m})$ be two ordered sets of $n\times n$ matrices. We say that $\mathrm{\mathbf{A}}_{[m]}$ and $\mathrm{\mathbf{B}}_{[m]}$ are simultaneously similar if there exists $\mathrm{\mathbf{X}}\in\mathrm{GL}(n)$ such that $\mathrm{\mathbf{A}}_{i}=\mathrm{\mathbf{X}}^{-1}\mathrm{\mathbf{B}}_{i}\mathrm{\mathbf{X}}$ for every $1\leq i\leq m$.

Let $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ be two $F$-PDs of the matrix multiplication tensor $\Phi_{m,p,n}$. For each $1\leq r\leq F$, define the matrices $\mathrm{\mathbf{M}}_{r}=\mathrm{\mathbf{W}}_{r}\mathrm{\mathbf{V}}_{r}\mathrm{\mathbf{U}}_{r}$ and $\mathrm{\mathbf{M}}^{\prime}_{r}=\mathrm{\mathbf{W}}^{\prime}_{r}\mathrm{\mathbf{V}}^{\prime}_{r}\mathrm{\mathbf{U}}^{\prime}_{r}$. If $\sigma((\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]}))$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ are (scaling+trace)-equivalent for some $\sigma\in\mathfrak{S}_{F}$, then from (7) we have

[TABLE]

In other words, $\sigma(\mathrm{\mathbf{M}}_{[F]})$ and ${\mathrm{\mathbf{M}}^{\prime}}_{[F]}$ are simultaneously similar.

We define a partial permutation of $\{1,\ldots,F\}$ as any injective function $\pi$ from $I\subseteq\{1,\ldots,F\}$ into $\{1,\ldots,F\}$. We say that $\pi$ coincides with the (total) permutation $\sigma\in\mathfrak{S}_{F}$ if $\pi(r)=\sigma(r)$ for every $r\in I$. If $\sigma$ is as in (12) and $\pi$ coincides with $\sigma$, then it is clear that

[TABLE]

The following notation will be useful for the description of the algorithm for computing $\sigma$. For $F^{\prime}\in\{0,\ldots,F\}$, we denote by $\mbox{Inj}(F^{\prime},F)$ the set of injective functions from $\{1,\ldots,F^{\prime}\}$ into $\{1,\ldots,F\}$. Each function of $\mbox{Inj}(F^{\prime},F)$ is seen as a subset of $\{1,\ldots,F^{\prime}\}\times\{1,\ldots,F\}$. The length of $\pi\in\mbox{Inj}(F^{\prime},F)$ is simply $\lvert\pi\rvert=F^{\prime}$, and the range of $\pi\in\mbox{Inj}(F^{\prime},F)$ is defined as $\mathrm{Range}(\pi)=\{\,\pi(r):1\leq r\leq F^{\prime}\,\}$.

The idea behind the algorithm to compute $\sigma$ is the following. First, we start from a partial permutation $\pi\in\mbox{Inj}(F^{\prime},F)$ with $F^{\prime}$ small. We check whether $\pi$ is susceptible to coincide with $\sigma$ by checking whether (13) is satisfied or not (see also Remark 4.4). If (13) is satisfied, then we try to extend $\pi$ to a larger partial permutation $\pi^{+}=\pi\cup\{(F^{\prime}+1,\ell)\}$ with $\ell\in\{1,\ldots,F\}\setminus\mathrm{Range}(\pi)$. We check again whether $\pi^{+}$ is susceptible to coincide with $\sigma$ according to (13). If this is the case, we repeat the process with $\pi^{+}$. Otherwise, we try other extensions $\pi\cup\{(F^{\prime}+1,\ell^{\prime})\}$. If all possible extensions $\pi\cup\{(F^{\prime}+1,\ell)\}$, $\ell\in\{1,\ldots,F\}\setminus\mathrm{Range}(\pi)$, have been tried and none of them coincides with $\sigma$, then we restart the process with the restriction $\pi^{-}=\pi|_{\{1,\ldots,F^{\prime}-1\}}\in\mbox{Inj}(F^{\prime}-1,F)$ and try to extend $\pi^{-}$ to $\pi^{-}\cup\{(F^{\prime},\ell)\}$ with $\ell\in\{1,\ldots,F\}\setminus\mathrm{Range}(\pi)$.

When we reach a full permutation $\pi\in\mbox{Inj}(F,F)=\mathfrak{S}_{F}$, then we can decide whether the permuted decomposition $\pi(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and the decomposition $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ are (scaling+trace)-equivalent using the procedure of the previous subsection. If they are, then we have found the correct permutation transformation between $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$. Otherwise, we continue to search for another permutation $\pi$.

When the algorithm terminates, if the two $F$-PDs are equivalent, the algorithm is guaranteed to give the corresponding scaling, trace and permutation transformations. If they are not equivalent, the algorithm will also detect it because all permutations $\pi\in\mathfrak{S}_{F}$ will be rejected: either because the partial permutation $\pi|_{\{1,\ldots,F^{\prime}-1\}}$ has been rejected previously in the algorithm, or because $\pi$ does not lead to (scaling+trace)-equivalent decompositions. Clearly, the computational savings (compared to trying all permutations) are interesting if most of the “incorrect” permutations $\pi$ are rejected in a early stage, i.e., $\pi|_{\{1,\ldots,F^{\prime}-1\}}$ is rejected for $F^{\prime}\ll F$. The computational aspects are discussed in the paragraphs below.

We have implemented the algorithm as the recursive function described in Algorithm 1. The recursive function must be called with (\pi,b)=\textnormal{{FnRecursive(}}\textnormal{\emph{\varnothing,\mbox{false}}}\textnormal{{)}}. If the output $b$ is true, then the two decompositions are equivalent and the permutation transformation is given by $\pi$. On the other hand, if $b$ is false, then the two $F$-PDs are not equivalent.

Remark 4.4*.*

Checking the simultaneous similarity of $\mathrm{\mathbf{A}}_{[m]}=(\mathrm{\mathbf{A}}_{1},\ldots,\mathrm{\mathbf{A}}_{m})$ and $\mathrm{\mathbf{B}}_{[m]}=(\mathrm{\mathbf{B}}_{1},\ldots,\mathrm{\mathbf{B}}_{m})$ can be approached by solving a linear system

[TABLE]

with unknown $\mathrm{\mathbf{X}}\in\mathbb{R}^{n\times n}$, and check whether there exists a solution $\mathrm{\mathbf{X}}$ that is invertible. However, this approach is not efficient and not robust to rounding errors. Therefore, we have used a different approach. Consider scalar coefficients $\alpha_{1},\ldots,\alpha_{m}\in\mathbb{R}$. A necessary condition for $\mathrm{\mathbf{A}}_{[m]}$ and $\mathrm{\mathbf{B}}_{[m]}$ to be simultaneously similar is that $\sum_{i=1}^{m}\alpha_{i}\mathrm{\mathbf{A}}_{i}$ and $\sum_{i=1}^{m}\alpha_{i}\mathrm{\mathbf{B}}_{i}$ have the same eigenvalues counted with multiplicity. By doing this for randomly generated sets of coefficients $\alpha_{1},\ldots,\alpha_{m}\in\mathbb{R}$, this gives a very efficient way to check the simultaneous similarity of $\mathrm{\mathbf{A}}_{[m]}$ and $\mathrm{\mathbf{B}}_{[m]}$ with high probability.  $\triangleleft$

Remark 4.5*.*

Strictly speaking, the use of Algorithm 1 supposes that Assumption 4.1 is satisfied. One could wonder whether we can still obtain some information from Algorithm 1 even if the assumption is not satisfied. The answer is yes. We modify the algorithm as follows. If condition $[\star]$ is satisfied, then instead of testing whether the $F$-PDs are (scaling+trace)-equivalent, we directly output $(\pi,\mbox{true})$ and exit the function. With this modified algorithm, if the call of the function (\pi,b)=\textnormal{{FnRecursive(}}\textnormal{\emph{\varnothing,\mbox{false}}}\textnormal{{)}} returns the value $b=\mbox{true}$, then we cannot say anything about the equivalence of the two $F$-PDs. However, if $b=\mbox{false}$, then we are sure that the two $F$-PDs are not equivalent.  $\triangleleft$

Numerical experiments for the algorithm described in this section are presented in Section 6.3. Regarding the complexity of the algorithm, the computation of the scaling and trace transformations relies only on solving linear systems of equations. The system (9) consists of $mpF$ equations with $(mp)^{2}+F$ variables. Because $F\geq mp$ (consequence of Theorem 4.4), the complexity of solving (9) is at most $\mathcal{O}([Fmp]^{3})$. Similarly, solving (11) requires at most $\mathcal{O}([F(pn+nm)]^{3})$. Therefore, the complexity of the (scaling+trace)-equivalence part of the algorithm is bounded by $\mathcal{O}([F\max\{mp,pn,nm\}]^{3})$.

The complexity of the permutation computation part is more difficult to evaluate. It is obviously bounded by $F!$. Hence, an upper bound for the global complexity of the algorithm is $\mathcal{O}(F![F\max\{mp,pn,nm\}]^{3})$. However, in all the numerical experiments we have performed (see Section 6.3), it appears that Algorithm 1 never reaches step $[\star]$ more than once. In fact, all the partial permutations $\pi$ for which the algorithm reaches step $[\spadesuit]$ satisfy $\lvert\pi\rvert\leq 9$ (see Table 2–Depth). In other words, whenever $\pi\in\mbox{Inj}(F^{\prime},F)$ could not lead to a correct permutation, then the algorithm detected it rapidly. Hence, in practice (for our numerical experiments), the computational complexity of the complete algorithm is $\mathcal{O}([F\max\{mp,pn,nm\}]^{3})$.

5 Characteristic polynomials and discretizable decompositions

Drawing upon the simultaneous similarity property (12) of equivalent decompositions, we introduce a simple necessary criterion for a decomposition of a matrix multiplication tensor to be equivalent to a discrete decomposition.

Definition 5.1**.**

A decomposition $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ is discretizable if it is equivalent to a discrete decomposition $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$. [Clearly, it is necessary and sufficient to require that $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ is only (scaling+trace)-equivalent to $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$.]

We refer the reader to Section 2 for the definition and relevance of discrete decompositions in the context of fast matrix multiplication. Numerical algorithms for computing polyadic decompositions of matrix multiplication tensors do not lead in general to solutions of this kind. The possibility to transform a general decomposition into a discrete one using invariance transformation opens the door to a new generation of algorithms to compute discrete solutions relying on a two-step approach (first compute a general decomposition and then discretize it). However, it is not clear when a decomposition can be discretized with invariance transformations so that the two-step approach may be inapplicable in some cases. The aim of this section is not to describe algorithms for transforming general decompositions into discrete decompositions but we propose a necessary criterion for a decomposition to be discretizable.

The criterion draws upon the observations made in Section 4.3: if $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ are (scaling+trace)-equivalent, then the families $\mathrm{\mathbf{M}}_{[F]}$ and ${\mathrm{\mathbf{M}}^{\prime}}_{[F]}$, defined by $\mathrm{\mathbf{M}}_{r}=\mathrm{\mathbf{W}}_{r}\mathrm{\mathbf{V}}_{r}\mathrm{\mathbf{U}}_{r}$ and $\mathrm{\mathbf{M}}^{\prime}_{r}=\mathrm{\mathbf{W}}^{\prime}_{r}\mathrm{\mathbf{V}}^{\prime}_{r}\mathrm{\mathbf{U}}^{\prime}_{r}$, are simultaneously similar (Definition 4.1). In particular, $\sum_{r=1}^{F}\beta_{r}\mathrm{\mathbf{M}}_{r}$ and $\sum_{r=1}^{F}\beta_{r}\mathrm{\mathbf{M}}^{\prime}_{r}$ are also similar for every coefficients $\beta_{r}\in\mathbb{R}$ (cf. Remark 4.4) and thus they have the same characteristic polynomial.

Assume that $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ is a discrete $F$-PD. Then $\mathrm{\mathbf{U}}^{\prime}_{r}\in(q\mathbb{Z})^{p\times m}$, $\mathrm{\mathbf{V}}^{\prime}_{r}\in(q\mathbb{Z})^{n\times p}$ and $\mathrm{\mathbf{W}}^{\prime}_{r}\in(q\mathbb{Z})^{m\times n}$ for some $q\in\mathbb{R}$. Hence, $\mathrm{\mathbf{M}}^{\prime}_{r}\in(q^{3}\mathbb{Z})^{m\times m}$ for every $1\leq r\leq F$. Let the coefficients $\beta_{r}$ in the paragraph above be integers. If we denote the characteristic polynomial of $\frac{1}{q^{3}}\sum_{r=1}^{F}\beta_{r}\mathrm{\mathbf{M}}_{r}$ by

[TABLE]

it is not hard to see that the coefficients $\alpha_{i}\in\mathbb{Z}$ for every $0\leq i<m$.

Definition 5.2**.**

Let the matrices $\mathrm{\mathbf{M}}_{r}$ be defined as above. We say that the decomposition $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ satisfies the discretizability criterion with parameter $q$ if for every integer coefficients $\beta_{r}$, $1\leq r\leq F$, the coefficients of the characteristic polynomial (14) satisfy $\alpha_{i}\in\mathbb{Z}$ for every $0\leq i<m$.

From the developments above, it is clear that satisfying the discretizability criterion with some parameter $q\in\mathbb{R}$ is a necessary condition for being discretizable. In the following section, we will see that most of the sample decompositions on which we have performed numerical experiments do not satisfy the discretizability criterion with $q=1$ or $q=1/2$ for tensors larger than the $2\times 2$ by $2\times 2$ case, contrasting with the abundance in the literature of discrete decompositions with $q=1$ or $q=1/2$ for these tensors (see also Section 1).

6 Numerical experiments

We have applied the results of Sections 4 and 5 on large sample sets of decompositions for matrix multiplication tensors up to the $m=p=n=3$ case. The goal is to get for the first time a view on the distributions of essentially unique decompositions and the distributions of discretizable decompositions: how many essentially unique decompositions do there exist? If two different decompositions are computed with a numerical algorithm, are they likely to be equivalent? Likely to be discretizable for some given $q$?

The way to obtain these samples is described in the next subsection. The reason we restrict to cases smaller than or equal to the $m=p=n=3$ case is explained in the next subsection as well. All computations were performed in Matlab. The computation-intensive part, namely, the generation of the samples, was executed on a Linux machine with 28 cores and 128 GBytes of RAM. The other computations were done on a laptop having 4 cores and 16 GBytes of RAM running Linux.

6.1 Computing polyadic decompositions

In the numerical experiments, we considered the six different cases $(m,p,n;F)$ summarized in Table 1 (first two columns), where $(m,p,n)$ is the size of the matrix multiplication tensor and $F$ is the number of rank-$1$ terms, i.e., we considered $F$-PDs of $\Phi_{m,p,n}$. For each $(m,p,n)$, the associated $F$ is the smallest $F$ for which we know there exists in the literature a decomposition of $\Phi_{m,p,n}$ with $F$ terms (see, e.g., [13, 25]).888Note that for the first four cases in Table 1, $F$ is equal to the rank of the associated tensor and thus cannot be decreased; see, e.g., [2, Chapter 15] for the $(1,2,1)$ and $(2,1,2)$ cases; for $(2,2,2)$, see [3, Theorem 11]; and for $(2,3,2)$, see [31]. For the $(3,2,3)$ case, the best known lower bound on the rank of $\Phi_{3,2,3}$ is $\mathrm{rank}(\Phi_{3,2,3})\geq 14$ (see, e.g., [3, Theorem 11]), and for $(3,3,3)$ the best known lower bound is $\mathrm{rank}(\Phi_{3,3,3})\geq 19$, shown by Bläser [8]. However, no $F$-term polyadic decompositions of $\Phi_{3,2,3}$ and $\Phi_{3,3,3}$ with $F<15$ and $F<23$ respectively are known for the moment. For each case, we want to obtain large sets of decompositions on which to apply the results of Sections 4 and 5.

Computing polyadic decompositions of matrix multiplication tensors is notoriously difficult (see, e.g., [25, 26] and references therein). Quite a few papers in the literature about tensor decompositions are devoted to this specific problem. For the numerical experiments of this paper, we have used the method proposed by Tichavský et al. [26] to compute $N_{s}=10\,000$ samples (decompositions) for the six cases listed in Table 1. For an alternative method, we refer the reader to [25]. See also [32]. We have used $\mathrm{\mathbf{U}}_{r,0}\in\mathbb{R}^{p\times m}$, $\mathrm{\mathbf{V}}_{r,0}\in\mathbb{R}^{n\times p}$ and $\mathrm{\mathbf{W}}_{r,0}\in\mathbb{R}^{m\times n}$ with entries chosen uniformly at random in $[-1,1]$ as initial iterates for Tichavský et al.’s method. The method does not always converge to a global minimum; hence we sometimes had to try more than one initial iterate to converge to an exact solution (the third and fourth columns give an idea of the effort required to compute the $N_{s}$ decompositions). In the end, we have at our disposal for each case $N_{s}$ samples of $F$-term polyadic decompositions of $\Phi_{m,p,n}$. We denote them by $({\mathrm{\mathbf{U}}^{\kappa}}_{[F]},{\mathrm{\mathbf{V}}^{\kappa}}_{[F]},{\mathrm{\mathbf{W}}^{\kappa}}_{[F]})$ with $\kappa\in\{1,\ldots,N_{s}\}$.

In the numerical computations, the tensors $\Phi\in\mathrm{Bil}(\mathbb{R}^{m\times p},\mathbb{R}^{p\times n};\mathbb{R}^{m\times n})$ are represented by the three-dimensional arrays $\tilde{\Phi}\in\mathbb{R}^{mp\times pn\times nm}$ obtained from the canonical identifications $\mathbb{R}^{m\times p}\cong\mathbb{R}^{mp}$, etc. Regarding floating-point arithmetic limitations, a sample $({\mathrm{\mathbf{U}}^{\kappa}}_{[F]},{\mathrm{\mathbf{V}}^{\kappa}}_{[F]},{\mathrm{\mathbf{W}}^{\kappa}}_{[F]})$ is considered as an $F$-PD of $\Phi_{m,p,n}$ if

[TABLE]

where $\Phi^{\kappa}$ is the tensor defined by (1) and (2) with $({\mathrm{\mathbf{U}}^{\kappa}}_{[F]},{\mathrm{\mathbf{V}}^{\kappa}}_{[F]},{\mathrm{\mathbf{W}}^{\kappa}}_{[F]})$.

Remark 6.1*.*

For matrix multiplication tensors larger than the $(3,3,3)$ case, it becomes very difficult to compute polyadic decompositions of these tensors: the global convergence of the algorithm decreases significantly while the cost for a single iteration of Tichavský et al.’s method grows as $\mathcal{O}([F(mp+pn+nm)]^{3})$. It becomes thus unrealistic to compute large sets of decompositions for these tensors.  $\triangleleft$

6.2 Discretizable decompositions

We start with the analysis of the discretizability property of the sample decompositions $({\mathrm{\mathbf{U}}^{\kappa}}_{[F]},{\mathrm{\mathbf{V}}^{\kappa}}_{[F]},{\mathrm{\mathbf{W}}^{\kappa}}_{[F]})$. For instance, we would like to find the decompositions that are not equivalent to a discrete decomposition with $q=1/2$ (Definition 2.1). To do this, we will apply the necessary criterion for discretizability with parameter $1/2$.

For every $\kappa\in\{1,\ldots,N_{s}\}$, let $p(t)=p(t;\beta_{1},\ldots,\beta_{F})$ be as in (14) where the $\beta_{r}$’s are randomly chosen integer coefficients. In our experiments, we use $16$ sets of coefficients sampled uniformly at random in $\{-5,-4,\ldots,5\}^{F}$, providing thus $16$ polynomials

[TABLE]

For each $\kappa\in\{1,\ldots,N_{s}\}$, we let

[TABLE]

where $\mathrm{round}(\alpha)$ is the closest integer to $\alpha$. The value of $\mbox{ND}_{\kappa}$ is thus a measure of how close are the polynomials $p^{\kappa}_{j}(t)$ to polynomials with integer coefficients. From the results of Section 5, a decomposition $({\mathrm{\mathbf{U}}^{\kappa}}_{[F]},{\mathrm{\mathbf{V}}^{\kappa}}_{[F]},{\mathrm{\mathbf{W}}^{\kappa}}_{[F]})$ for which $\mbox{ND}_{\kappa}$ is (significantly) nonzero is not equivalent to a discrete decomposition with $q=1/2$.

The histograms in Figure 2 show the distribution of decompositions based on the value of $\mbox{ND}_{\kappa}$. More precisely, each bar of the histograms represents the number of decompositions $({\mathrm{\mathbf{U}}^{\kappa}}_{[F]},{\mathrm{\mathbf{V}}^{\kappa}}_{[F]},{\mathrm{\mathbf{W}}^{\kappa}}_{[F]})$ with $\mbox{ND}_{\kappa}$ in the corresponding range. For the $(1,2,1)$, and $(2,2,2)$ cases, we observe that, for all the decompositions, the polynomials $p^{\kappa}_{j}(t)$ have integer coefficients (within a very small tolerance). Hence, $100\%$ of the decompositions satisfy the discretizability criterion with $q=1/2$. This is not surprising since all the $2$-PDs (resp. $7$-PDs) of $\Phi_{1,2,1}$ (resp. $\Phi_{2,2,2}$) are equivalent (see [4] and Remark 6.2), and $\Phi_{1,2,1}$ (resp. $\Phi_{2,2,2}$) admits a discrete decomposition with $q=1$.999For $\Phi_{1,2,1}$, take, e.g.,

$\displaystyle\mathrm{\mathbf{U}}_{1}=[1,0],\quad\mathrm{\mathbf{V}}_{1}=[1,0]^{\top},\quad\mathrm{\mathbf{W}}_{1}=1,$

$\displaystyle\mathrm{\mathbf{U}}_{2}=[0,1],\quad\mathrm{\mathbf{V}}_{2}=[0,1]^{\top},\quad\mathrm{\mathbf{W}}_{2}=1,$

(see also Remark 6.2). For $\Phi_{2,2,2}$, take, e.g., Strassen’s algorithm.

In contrast, for the $(2,1,2)$, $(2,3,2)$, $(3,2,3)$ and $(3,3,3)$ cases, we observe that most of the decompositions do not satisfy the necessary criterion for discretizability with $q=1/2$. This implies that most of the decompositions are not equivalent to a discrete decompositions with parameter $q=1/2$. This last observation has to be put in contrast with the abundance of decompositions $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ for which the entries of $\mathrm{\mathbf{U}}_{r}$, $\mathrm{\mathbf{V}}_{r}$ and $\mathrm{\mathbf{W}}_{r}$ belong to $\{0,\pm 1/2,\pm 1\}$ in the literature [13, 5, 26, 24].

Further experiments can be conducted to investigate the discretizability of the decompositions with respect to other parameters $q$. However, due to space limitations, we do not present them in this paper.

6.3 Equivalence classes of decompositions

In the previous subsection, we have seen that, except for the $(1,2,1)$ and $(2,2,2)$ cases, most of the decompositions are not equivalent to a discrete decomposition with coefficients in $\{0,\pm 1/2,\pm 1\}$. In this section, we will analyze the pairwise equivalence of the decompositions. This will reveal the distributions of the equivalence classes among the sample sets of decompositions. Therefore, we use the algorithm developed in Section 4.

First, in order to apply Algorithm 1, we need to ensure that Assumption 4.1 is satisfied for every decomposition $({\mathrm{\mathbf{U}}^{\kappa}}_{[F]},{\mathrm{\mathbf{V}}^{\kappa}}_{[F]},{\mathrm{\mathbf{W}}^{\kappa}}_{[F]})$, $\kappa\in\{1,\ldots,N_{s}\}$. Therefore, for each decomposition $({\mathrm{\mathbf{U}}^{\kappa}}_{[F]},{\mathrm{\mathbf{V}}^{\kappa}}_{[F]},{\mathrm{\mathbf{W}}^{\kappa}}_{[F]})$, we have computed (using Theorem 4.3) the clustering vector of the decomposition defined as the vector $[\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{U}}}^{\kappa}),\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{V}}}^{\kappa}),\mathrm{cl}_{\oplus}(\tilde{\mathrm{\mathbf{W}}}^{\kappa})]$. The results are summarized in Figure 3. As we can see, for each case, $100\%$ of the decompositions have at least one matrix $\tilde{\mathrm{\mathbf{U}}}^{\kappa}$, $\tilde{\mathrm{\mathbf{V}}}^{\kappa}$ or $\tilde{\mathrm{\mathbf{W}}}^{\kappa}$ with clustering number equal to one and thus satisfy Assumption 4.1.

We can thus apply Algorithm 1 to check the equivalence between pairs of decompositions for the different cases. The results are gathered in Table 2. We observe that for the $(1,2,1)$ and $(2,2,2)$ cases, every decompositions are pairwise equivalent (see also Remark 6.2). As a consequence, it is not surprising that all the decompositions are discretizable. This situation is more surprising for $(2,1,2)$, $(2,3,2)$, $(3,2,3)$ and $(3,3,3)$ cases. For these cases, the decompositions seem to be pairwise equivalent with probability zero.

The third column of Table 2 gives the average computation time to check the equivalence between two decompositions. We observe that the algorithm takes no more than $30~{}\text{ms}$. In comparison, for the $(3,3,3)$ case for example, the naive method (testing all possible permutations) would have required to test condition $[\clubsuit]$ in Algorithm 1, which has a complexity of $\mathcal{O}([F\max\{mp,pn,nm\}]^{3})$, $23!=2.59\cdot 10^{22}$ times.

Remark 6.2*.*

It can be shown that all the decompositions of $\Phi_{1,2,1}$ and $\Phi_{2,2,2}$ (respectively) are pairwise equivalent, which corroborates the results of the numerical experiments using Algorithm 1. For the $(2,2,2)$ case, we refer the reader to [4]. For the $(1,2,1)$ case, let $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ and $({\mathrm{\mathbf{U}}^{\prime}}_{[F]},{\mathrm{\mathbf{V}}^{\prime}}_{[F]},{\mathrm{\mathbf{W}}^{\prime}}_{[F]})$ be two $2$-PDs of $\Phi_{1,2,1}$. Observe that $\Phi_{1,2,1}$ maps $2$-dimensional vectors to their scalar product and thus can be represented with the identity matrix:

[TABLE]

Since $(\mathrm{\mathbf{U}}_{[F]},\mathrm{\mathbf{V}}_{[F]},\mathrm{\mathbf{W}}_{[F]})$ is a decomposition, it is not hard to see that

[TABLE]

where we remind that $\tilde{\mathrm{\mathbf{U}}}$ and $\tilde{\mathrm{\mathbf{V}}}$ are defined as (6). Using a scaling transformation, we may assume that $\mathrm{\mathbf{W}}_{1}=\mathrm{\mathbf{W}}_{2}=1$. Hence $\tilde{\mathrm{\mathbf{U}}}$ and $\tilde{\mathrm{\mathbf{V}}}^{\top}$ are inverses of each other, and so are $\tilde{\mathrm{\mathbf{U}}}^{\prime}$ and $\tilde{\mathrm{\mathbf{V}}}^{\prime\top}$. Then let $\mathrm{\mathbf{P}}=\tilde{\mathrm{\mathbf{U}}}\tilde{\mathrm{\mathbf{V}}}^{\prime\top}$ and observe that

[TABLE]

Hence, we have found a trace transformation with $\mathrm{\mathbf{Q}}=\mathrm{\mathbf{R}}=1\in\mathrm{GL}(1)$ between the two decompositions.  $\triangleleft$

7 Conclusions

In this paper, we have described an algorithm for efficiently deciding whether two decompositions of a given matrix multiplication tensor are equivalent through invariance transformations. We have introduced the notion of clustering number of a matrix and we have demonstrated the correctness of the algorithm provided some conditions on the clustering number of the factor matrices of the decompositions are satisfied. This condition was satisfied for $100\%$ of the numerical samples on which we have applied our algorithm.

The analysis of the equivalence classes of decompositions is relevant in the context of fast matrix multiplication as it sheds light on the diversity of essentially unique fast matrix multiplication algorithms. In the numerical experiments we have performed, it appears that two decompositions are equivalent with probability zero (except for the multiplication of $1\times 2$ by $2\times 2$ matrices and the multiplication of $2\times 2$ matrices for which we can prove the essential uniqueness of their decompositions) indicating that there are many essentially different algorithms for the fast multiplication of $3\times 3$ matrices for example.

Drawing upon the observation that decompositions with coefficients in a discrete set provide fast matrix multiplication with better performance, we have also provided a necessary criterion for a decomposition to be equivalent to a decomposition with these properties. We have applied the criterion on numerical samples and observed that the majority of the decompositions do not satisfy the criterion for being equivalent to a decomposition with coefficients in, e.g., $\{0,\pm 1\}$ or $\{0,\pm 1/2,\pm 1\}$.

Acknowledgments

The authors would like to thank Samuel Fiorini and Gwenaël Joret for insightful discussions on the link between the clustering number and matroids.