Jean-Guillaume Dumas
email: jean-guillaume.dumas@univ-grenoble-alpes.fr

Clément Pernet
email: clement.pernet@univ-grenoble-alpes.fr

Alexandre Sedoglavic
email: alexandre.sedoglavic@univ-lille.fr

Some fast algorithms multiplying a matrix by its adjoint

We present a non-commutative algorithm for the multiplication of a 2 × 2 block-matrix by its adjoint, defined by a matrix ring anti-homomorphism. This algorithm uses 5 block products (3 recursive calls and 2 general products)over C or in positive characteristic. The resulting algorithm for arbitrary dimensions is a reduction of multiplication of a matrix by its adjoint to general matrix product, improving by a constant factor previously known reductions. We prove also that there is no algorithm derived from bilinear forms using only four products and the adjoint of one of them. Second we give novel dedicated algorithms for the complex field and the quaternions to alternatively compute the multiplication taking advantage of the structure of the matrix-polynomial arithmetic involved. We then analyze the respective ranges of predominance of the two strategies. Finally we propose schedules with low memory footprint that support a fast and memory efficient practical implementation over a prime field.

Introduction

Volker Strassen's algorithm [START_REF] Strassen | Gaussian elimination is not optimal[END_REF], with 7 recursive multiplications and 18 additions, was the first sub-cubic time algorithm for matrix product, with a cost of O n 2.81 . Summarizing the many improvements which have happened since then, the cost of multiplying two arbitrary n × n matrices over a ring R will be denoted by MM R ω (n) = O(n ω ) ring operations where 2 < ω ≤ 3 is any feasible exponent for this operation (see [START_REF] Gall | Powers of tensors and fast matrix multiplication[END_REF] for the best theoretical estimates of ω known to date).

We consider here the computation of the product of a matrix by transpose A • A or by its conjugate transpose A • A , which we handle in a unified way as the product A • φ (A) where φ is a matrix anti-homomorphism. In the rest of the paper, φ (A) will be referred to as the adjoint of A. For this computation, the natural divide and conquer algorithm, splitting the matrices in four quadrants, will use 6 block multiplications (as any of the two off-diagonal blocks can be recovered from the other one).We propose instead a new algorithm using only 5 block multiplications, for any antihomomorphism φ, provided that the base ring supports the existence of skew unitary matrices.

2

For this product, the best previously known cost bound was equivalent to 2 2 ω -4 MM ω (n) over any field (see [6, § 6.3.1]). With our algorithm, this product can be computed in a cost equivalent to 2 2 ω -3 MM ω (n) ring operations when there exists a skew-unitary matrix. Our algorithm is derived from the class of Strassen-like algorithms multiplying 2 × 2 matrices in 7 multiplications. Yet it is a reduction of multiplying a matrix by its transpose to general matrix multiplication, thus supporting any admissible value for ω. By exploiting the symmetry of the problem, it requires about half of the arithmetic cost of general matrix multiplication when ω is log 2 7.

This paper extends the results of [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF] with the following improvements:

1. we generalize the case of the transposition in [7, Algorithm 2] to arbitrary antihomomorphism, including the Hermitian transposition. 2. Our algorithm uses 5 multiplications and the (hermitian) transpose of one these blocks. In [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF] a Gröbner basis parameterization is used to search for algorithms, or prove by exhaustive search that there are no better algorithm, in the Strassen orbit. We partially address here the more general result that there is no algorithm derived from bilinear forms, with fewer products, by proving the inexistence of an algorithm with four products and the (hermitian) transpose of one of them. 3. In [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF] the algorithm is shown to be efficient over C, for a range of matrix multiplication exponents (including all the feasible ones), and for any positive characteristic field, unconditionally. We extend this analysis to the case for the Hermitian transpose: while our five-products algorithm is unusable due to the inexistence of skew unitary matrices over C, we propose a 2M algorithm, adapted from the 3M algorithm for the product of complex matrices. [START_REF] Brillhart | Note on representing a prime as a sum of two squares[END_REF]. Finally, we propose novel dedicated algorithms for the multiplication of a matrix by its transpose or conjugate transpose over the algebra of quaternions (over R or any commutative field), improving on the dominant term of the state of the art complexity bounds for these problems. After a introducing the terminology in Section 2, we will present in Section 3 the main recursive algorithm computing the product of a matrix by its adjoint in 5 block products provided that a skew unitary matrix is given. We survey in Section 4 the most classical instances for the base field to support the existence of skew unitary matrices. We then investigate in Section 5 the minimality of five products for the computing the product of a 2 × 2 matrix by its hermitian transpose: applying de Groote's technique enables us to state this result partially, for all algorithms using up to one symmetry between a product and its adjoint. Section 6 explores alternative approaches offered by the structure of polynomial arithmetic, when the field is an extension. This includes a new 2M algorithm in Section 6.1 and new algorithms over the algebra of quaternions in Section 6.2. Lastly, we discuss on an implementation of the recursive algorithm for the product of a matrix by its transpose in Section 7.

Preliminaries

To unify the notion of transposition and conjugate transposition, we use the formalism of antihomomorphisms and of involutive antihomomorphisms as recalled in the following definitions. Definition 1. Let R, S be two rings, φ : R → S is a ring antihomomorphism if and only if, for all (x, y) in R × S:

φ(1 R ) = 1 S , (1a)
φ(x + y) = φ(x) + φ(y), (1b)
φ(xy) = φ(y)φ(x). (1c)
From this, one can define a matrix antihomomorphism by induction, as shown in Definition 2.

Definition 2. Over a ring R, an involutive matrix antihomomorphism is a family of applications φ m,n : R m×n → R n×m for all (m, n) in N 2 satisfying for additional ( , k) in N 2 and for all A and A in R m×n , for all M in R n×k , for all B in R m×k , for all C in R ×n and for all D in R ×k the following relations:

φ m,n • φ n,m = I, (2a)
φ m,n (A + A ) = φ m,n (A) + φ m,n (A ), (2b)
φ m,k (A • M ) = φ n,k (M ) • φ m,n (A), (2c)
φ m+ ,n+k A B C D = φ m,n (A) φ ,n (C) φ m,k (B) φ ,k (D) . ( 2d
)
For the convenience, we will denote all applications of this family by φ, as the dimensions are clear from the context. This definition implies the following: Lemma 3. For all A in R m×n let B, be φ (A). Then for all suitable (i, j) the coefficient b ij is φ(a ji ).

Proof. By induction, using Equation (2d): if

m = n = 1, then φ (A) = [φ(a 11 )].
Then assume the property is true for all A ∈ R m×n with m, n ≤ N , and consider a matrix A in R (N +1)×(N +1) . Applying Equation (2d) on the block decomposition A = A 11 a 12 a 21 a 22 where A 11 is in R N ×N yields the relations:

φ (a) = φ (A 11 ) φ (a 21 ) φ (a 12 ) φ(a 22 ) = [φ(a ji )] ij (3)
by induction hypothesis. The case of matrices in R m×(N +1) and R (N +1)×n is dealt with similarly, using 0-dimensional blocks a 21 or a 12 respectively.

Lemma 4. For all α in R and for all A in R m×n , φ (αA) = φ (A)φ (α).

Proof. By Lemma 3, φ (αI m ) = φ(α)I m = I m φ(α). Then by Equation (2c), the relations φ (αA

) = φ ((αI m )A) = φ (A)φ (αI m ) = φ (A)φ(α) hold.
The following Lemma 5 shows that Definition 2 is a natural extension of a ring antiendomorphism for matrices. Lemma 5. An involutive matrix antihomomorphism is a ring antiendomorphism on its base ring (seen as the ring of 1×1 matrices).

Proof. Equations (2b) and (2c) directly imply Equations (1b) and (1c) respectively when m = n = k = 1. Then, we have that φ(1) = φ(1) • 1. Therefore φ(φ(1)) = φ(φ(1) • 1) and 1 = φ(1)φ(φ(1)) = φ(1) • 1 by Equations (2a) and (2c). This right hand side is equal to that of the first equation, thus proving the equality of the left hand sides and Equation (1a).

Definition 2 gathers actually all the requirements for our algorithm to work in classical hermitian or non-hermitian cases: Examples 6. For matrices over a commutative ring,

• the matrix transpose with φ (A) = A and

• the matrix conjugate transpose, φ (A) = A H ,
are two examples of matrix anti-homomorphisms. However, for instance, transposition over the quaternions is a counter-example as the non-commutativity implies there that in general

(A • B) = B • A .
• L ij = a ij for i ≥ j and L ij = 0 otherwise, • U ij = a ij for i ≤ j and U ij = 0 otherwise. Lemma 9. If φ (A) = A in R n×n , then Up (A) = φ (Low (A)). Proof. Applying Lemma 3, the coefficients u ij of U = φ (Low (A)) for 0 < i ≤ j satisfy u ij = φ(a ji ). Now if φ (A) = A, we have u ij = a ij for 0 < i ≤ j and u ij = 0 otherwise, as φ (0) = 0, by Equation (2b). Hence U = Up (A).
Definition 10 (Skew-unitary). A matrix Y in R n×n is skew-unitary relatively to a matrix antihomomorphism φ if the following relation holds:

Y • φ (Y ) = -I n . (4)
For the cost analysis, we will also need the following variant of the Master Theorem, reflecting the constant in the leading term of the computed cost bound.

Lemma 11. Let T (n) be defined by the recurrence

T (n) = aT (n/2) + b n 2 α + o(n α ), where 0 ≤ log 2 a < α. Then T (n) = b 2 α -a n α + o(n α ). Proof. T (n) = a log 2 n T (1) + log 2 (n)-1 i=0 a i b n 2 i+1 α + o n 2 i α = n log 2 a T (1) + b 2 α n α log 2 (n)-1 i=0 a 2 α i + o(n α ) = b 2 α -a n α + o(n α ).
3 An algorithm for the product of a matrix by its adjoint with five multiplications

We now show how to compute the product of a matrix by its adjoint with respect to an involutive antihomomorphism in only 5 recursive multiplications and 2 multiplications by any skew-unitary matrix. This is a generalization of [7, Algorithm 2] for any involutive antihomomorphism. We next give Algorithm 12 for even dimensions. In case of odd dimensions, padding or static/dynamic peeling can always be used [START_REF] Boyer | Memory efficient scheduling of Strassen-Winograd's matrix multiplication algorithm[END_REF].

Theorem 13. Algorithm 12 is correct. Moreover, if any two n × n matrices over a ring R can be multiplied in MM R ω (n) = O(n ω ) ring operations for ω > 2, and if there exist a skew-unitary matrix which can be multiplied to any other matrix in o(n ω ) ring operations then Algorithm 12 requires fewer than 2 2 ω -3 MM R ω (n) + o(n ω ) ring operations. Proof. For the cost analysis, Algorithm 12 is applied recursively to compute three products P 1 , P 2 and P 7 , while P 4 and P 5 are computed in MM R ω (n) using the general matrix multiplication algorithm. The second hypothesis is that applying the skew-unitary matrix Y to a n × n matrix costs Y (n) = o(n ω ). Then applying Remark 15 thereafter, the cost T (n) of Algorithm 12 satisfies:

T (n) ≤ 3T (n/2) + 2MM R ω (n/2) + 2Y (n) + (7.5)(n/2) 2 + o n 2 (5)
and T (4) is a constant. Thus, by Lemma 11:

T (n) ≤ 2 2 ω -3 MM R ω (n) + o(n ω ). ( 6
)
Now for the correction, by Equation (2d), we have to show that the result of Algorithm 12 is indeed:

Low (A • φ (A)) = Low A • φ(A11) φ(A21) φ(A12) φ(A22) = Low(A11•φ(A11)+A12•φ(A12)) × A21•φ(A11)+A22•φ(A12) Low(A21•φ(A21)+A22•φ(A22))
Algorithm 12 Product of a matrix by its adjoint Input: A ∈ R m×n (with even m and n for the sake of simplicity); Input: φ an involutive matrix antihomomorphism;

Input: Y ∈ R n 2 × n 2 skew-unitary for φ. Output: Low (A • φ (A)). Split A = A11 A12 A21 A22 where A 11 is in R m 2 × n 2
4 additions and 2 multiplications by Y :

S 1 ← (A 21 -A 11 ) • Y S 2 ← A 22 -A 21 • Y S 3 ← S 1 -A 22 S 4 ← S 3 + A 12
3 recursive (P1, P2, P5) and 2 general products (P3, P4):

Low (P 1 ) ← Low (A 11 • φ (A 11 )) Low (P 2 ) ← Low (A 12 • φ (A 12 )) P 3 ← A 22 • φ (S 4 ) P 4 ← S 1 • φ (S 2 ) Low (P 5 ) ← Low (S 3 • φ (S 3 ))
3 half additions and 2 complete additions:

Low (U 1 ) ← Low (P 1 )+Low (P 5 ) Low (U 3 ) ← Low (P 1 )+Low (P 2 ) Up(U 1 ) ← φ (Low (U 1 ))

Forms the full matrix U1

U 2 ← U 1 + P 4 , U 4 ← U 2 + P 3 , Low (U 5 ) ← Low (U 2 ) + Low (φ (P 4 )). return Low(U3) U4
Low(U5) .

First, we have that:

Low (U 3 ) = Low (P 1 ) + Low (P 2 ) = Low (A 11 • φ (A 11 ) + A 12 • φ (A 12 )). (7)
Second, as Y is skew-unitary, then we have that Y • φ (Y ) = -I n 2 . Also, by Equations (2b) and (2c), φ (S 2 ) = φ (A 22 ) -φ (Y ) • φ (A 21 ). Then, denote by R 1 the product:

R 1 = A 11 • Y • φ (S 2 ) = A 11 • Y • (φ (A 22 ) -φ (Y ) • φ (A 21 )) = A 11 • (Y • φ (A 22 ) + φ (A 21 )). (8)
Further, by Equations (2a) and (2c), we have that

P 1 = A 11 • φ (A 11 ), P 2 = A 12 • φ (A 12
), and P 5 = S 3 • φ (S 3 ) are invariant under the action of φ. So are therefore, U 1 = P 1 + P 5 , U 3 = P 1 + P 2 and U 5 = U 1 + (P 4 + φ (P 4 )). By Lemma 9, it suffices to compute Low (U 1 ) and, if needed, we also have Up (U 1 ) = φ (Low (U 1 )).

Then, as

S 3 = S 1 -A 22 = (A 21 -A 11 ) • Y -A 22 = -S 2 -A 11 • Y , and
φ (φ (S 2 )) = S 2 by Equation (2a), we have that:

U 1 = P 1 + P 5 = A 11 • φ (A 11 ) + S 3 • φ (S 3 ) = A 11 • φ (A 11 ) + (S 2 + A 11 • Y ) • (φ (S 2 ) + φ (Y ) • φ (A 11 )) = S 2 • φ (S 2 ) + φ (R 1 ) + R 1 . (9)
Also, denote

R 2 = A 21 • Y • φ (A 22
), so that:

S 2 • φ (S 2 ) = (A 22 -A 21 • Y ) • (φ (A 22 ) -φ (Y ) • φ (A 21 )) = A 22 • φ (A 22 ) -A 21 • φ (A 21 ) -R 2 -φ (R 2 ). ( 10
)
Furthermore, from Equation (8):

R 1 + P 4 = R 1 + S 1 • φ (S 2 ) = R 1 + (A 21 -A 11 ) • Y • (φ (A 22 ) -φ (Y ) • φ (A 21 )) = A 21 • Y • φ (A 22 ) + A 21 • φ (A 21 ) = R 2 + A 21 • φ (A 21 ). (11)
This shows, from Equations ( 9) to [START_REF] Friedrich Groote | On the complexity of quaternion multiplication[END_REF], that:

U 5 = U 1 + P 4 + φ (P 4 ) = S 2 • φ (S 2 ) + φ (R 1 ) + R 1 + P 4 + φ (P 4 ) = A 22 • φ (A 22 ) + (-1 + 2)A 21 • φ (A 21 ). (12)
Third, the last coefficient U 4 of the result is obtained from Equations ( 11) and [START_REF] Friedrich Groote | On varieties of optimal algorithms for the computation of bilinear mappings I. The isotropy group of a bilinear mapping[END_REF]:

U 4 = U 2 + P 3 = U 1 + P 4 + P 3 = A 22 • φ (A 22 ) -A 21 • φ (A 21 ) -R 2 -φ (R 2 ) + R 1 + φ (R 1 ) + P 4 + P 3 = A 22 • φ (A 22 ) -φ (R 2 ) + φ (R 1 ) + P 3 (13)
since by Equation [START_REF] Friedrich Groote | On the complexity of quaternion multiplication[END_REF]

, R 1 + P 4 = R 2 + A 21 • φ (A 21 ). Now P 3 = A 22 • φ (S 4 ) = A 22 • φ ((A 21 -A 11 ) • Y + A 12 -A 22 ) = φ (R 2 ) -φ (R 1 ) + A 21 • φ (A 11 ), (14)
Hence

U 4 = A 22 • φ (A 12 ) + A 21 • φ (A 11 )
To our knowledge, the best previously known result was with a 2 2 ω -4 factor instead, see e.g. [6, § 6.3.1]. Table 1 summarizes the arithmetic complexity bound improvements.

Examples 14. In many cases, applying the skew-unitary matrix Y to a n × n matrix costs only yn 2 for some constant y depending on the base ring. If the ring is the complex field C or satisfies the conditions of Proposition 16, there is a square root i of -1.

Setting Y = i I n/2 yields Y (n) = n 2 . Otherwise, we Problem Alg. O n 3 O n log 2 (7) O(n ω ) A • φ (A) ∈ F n×n [6] n 3 2 3 MM log 2 (7) (n) 2 2 ω -4 MM ω (n) Alg. 12 0.8n 3 1 2 MM log 2 (7) (n) 2 2 ω -3 MM ω (n)
Table 1: Arithmetic complexity bounds leading terms.

show in Section 4 that in characteristic p ≡ 3 mod 4, Proposition 17 produces Y equal to a b -b a ⊗ I n/2 for which Y (n) = 3n 2 . As a sub-case, the latter can be improved when p ≡ 3 mod 8: then, Lemma 18 shows that -2 is a square. Therefore, in this case set a = 1 and b ≡ √ -2 mod p such that one multiplication is saved. Then the relation

a 2 + b 2 = -1 there yields Y = 1 √ -2 - √ -2 1 ⊗ I n/2 for which Y (n) = 2n 2 .
Remark 15. Each recursive level of Algorithm 12 is composed of 9 block additions. An exhaustive search on all symmetric algorithms in the orbit of that of Strassen (via a Gröbner basis parameterization [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF]) showed that this number is minimal in this class of algorithms. Note also that 3 out of these 9 additions in Algorithm 12 involve symmetric matrices and are therefore only performed on the lower triangular part of the matrix. Overall, the number of scalar additions is 6n 2 + 3/2n(n + 1) = 15/2n 2 + 1.5n, nearly half of the optimal in the non-symmetric case [5, Theorem 1].

To further reduce the number of additions, a promising approach is that undertaken in [START_REF] Karstadt | Matrix multiplication, a little faster[END_REF][START_REF] Beniamini | Faster matrix multiplication via sparse decomposition[END_REF]. This is however not clear to us how to adapt our strategy to their recursive transformation of basis.

Rings with skew unitary matrices

Algorithm 12 requires a skew-unitary matrix. Unfortunately there are no skewunitary matrices over R, nor Q for φ the transposition, nor over C for φ the Hermitian transposition (there -1 cannot be a sum of real squares for a diagonal element of Y • φ (Y )). Hence, Algorithm 12 provides no improvement in these cases. In other domains, the simplest skew-unitary matrices just use a square root of -1 while others require a sum of squares.

Over the complex field

Algorithm 12 is thus directly usable over C n×n with φ (A) = A and

Y = i I n 2 in C n 2 × n 2 .
When complex numbers are represented in Cartesian form, as a pair of real numbers, the multiplications by Y = i I n 2 are essentially free since they just exchange the real and imaginary parts, with one sign flip.

As mentioned, for the conjugate transposition, φ (A) = A H , on the contrary, there are no candidate skew-unitary matrices and we for now report no im-provement in this case using this approach (but another one does as shown in Section 6.1). Now, even though over the complex the product of a matrix by its conjugate transpose is more widely used, there are some applications for the product of a matrix by its transpose, see for instance [START_REF] Baboulin | A parallel distributed solver for large dense symmetric systems: Applications to geodesy and electromagnetism problems[END_REF]. This is reflected in the blas api, where both routines zherk and zsyrk are offered.

Rings where negative one is a square and φ (A) = A

Over some rings , square roots of -1 can also be elements of the base field, denoted i in R again. There, Algorithm 12 only requires some pre-multiplications by this square root (with also

Y = i I n 2 ∈ R n 2 × n
2 ), but within the ring. Further, when the ring is a field in positive characteristic, the existence of a square root of minus one can be characterized, as shown in Proposition 16, thereafter.

Proposition 16. Fields with characteristic two, p satisfying p ≡ 1 mod 4, or finite fields that are an even extension, contain a square root of -1. = ±1 and then the square of the latter must be -1. If the finite field F is of cardinality p 2k , then, similarly, there exists elements

Proof. If p = 2, then 1 = 1 2 = -1. If p ≡ 1 mod 4,
x p k -1 2 p k +1 2
different from ±1 and then the square of the latter must be -1.

Any field with positive characteristic and φ (A) = A

Actually, we show that Algorithm 12 can also be run without any field extension, even when -1 is not a square: form the skew-unitary matrices constructed in Proposition 17, thereafter, and use them directly as long as the dimension of Y is even. Whenever this dimension is odd, it is always possible to pad with zeroes so that A

• A = ( A 0 ) • A 0 . Proposition 17. Let F p k be a field of characteristic p, there exists (a, b) in F 2 p such that the matrix: a b -b a ⊗ I n = a In b In -b In a In in F 2n×2n p (15)
is skew-unitary for the transposition.

Proof. Using the relation

a In b In -b In a In a In b In -b In a In = (a 2 + b 2 ) I 2n , (16)
it suffices to prove that there exist a, b such that

a 2 + b 2 = -1. In character- istic 2, a = 1, b = 0 is a solution as 1 2 + 0 2 = -1.
In odd characteristic, there are p+1 2 distinct square elements x i 2 in the base prime field. Therefore, there are p+1 2 distinct elements -1 -x i 2 . But there are only p distinct elements in the base field, thus there exists a couple (i, j) such that -1 -

x i 2 = x j 2 [18, Lemma 6].
To further improve the running time of multiplications by a skew-unitary matrix in this case, one could set one of the squares to be 1. This is possible if -2 is a square, for instance when p ≡ 3 mod 8:

Lemma 18. If p ≡ 3 mod 8 then -2 is a square modulo p. Proof. Using Legendre symbol, -2 p = -1 p 2 p = (-1) p-1 2 (-1) p 2 -1 8 = (-1)(-1) = 1
Now, Proposition 17 shows that skew-unitary matrices do exist for any field with positive characteristic. For Algorithm 12, we need to build them mostly for p ≡ 3 mod 4 (otherwise use Proposition 16).

For this, without the extended Riemann hypothesis (erh), it is possible to use the decomposition of primes into squares:

1. Compute by enumeration a prime r = 4pk + (3 -1)p -1, so that both relations r ≡ 1 mod 4 and r ≡ -1 mod p hold;

2. Thus, the methods of [START_REF] Brillhart | Note on representing a prime as a sum of two squares[END_REF] allow one to decompose any prime into squares and give a couple

(a, b) in Z 2 such that a 2 + b 2 = r. Finally, this gives a 2 + b 2 ≡ -1 mod p.
By the prime number theorem the first step is polynomial in log(p), as is the second step (square root modulo a prime, denoted sqrt, has a cost close to exponentiation and then the rest of Brillhart's algorithm is gcd-like). In practice, though, it is faster to use the following Algorithm 19, even though the latter has a better asymptotic complexity bound only if the erh is true.

Algorithm 19 SoS: Sum of squares decomposition over a finite field

Input: p ∈ P\{2}, k ∈ Z. Output: (a, b) ∈ Z 2 , s.t. a 2 + b 2 ≡ k mod p. 1: if k p = 1 then k is a square mod p 2:
return (sqrt(k), 0).

3: else

Find smallest quadratic non-residue

4: s ← 2; while s p == 1 do s ← s + 1 5: c ← sqrt(s -1)
s -1 must be a square Proof. If k is square then the square of one of its square roots added to the square of zero is a solution. Otherwise, the lowest quadratic non-residue (lqnr) modulo p is one plus a square b 2 (1 is always a square so the lqnr is larger than 2). For any generator of Z p , quadratic non-residues, as well as their inverses (s is invertible as it is non-zero and p is prime), have an odd discrete logarithm. Therefore the multiplication of k and the inverse of the lqnr must be a square a 2 . This means that the relation k = a 2 1 + b 2 = a 2 + (ab)

6: r ← ks -1 mod p 7: a ← sqrt(r) Now k ≡ a 2 s ≡ a 2 (1 + c 2 ) mod
2 holds. Now for the running time, under the erh, [START_REF] Wedeniwski | Primality tests on commutator curves[END_REF]Theorem 6.35] shows that the lqnr should be lower than 3 log 2 (p)/2 -44 log(p)/5 + 13. From this, the expected number of Legendre symbol computations is O log 2 (p) and this dominates the modular square root computations.

Remark 21. Another possibility is to use randomization: instead of using the lowest quadratic non-residue ( lqnr), randomly select a non-residue s, and then decrement it until s -1 is a quadratic residue (1 is a square so this will terminate). In practice, the running time seems very close to that of Algorithm 19 anyway, see, e.g. the implementation in Givaro rev. 7bdefe6, https: // github. com/ linbox-team/ givaro . Also, when computing t sum of squares modulo the same prime, one can compute the lqnr only once to get all the sum of squares with an expected cost bounded by O log 3 (p) + tlog 2 (p) .

Remark 22. Except in characteristic 2 or in algebraic closures, where every element is a square anyway, Algorithm 19 is easily extended over any finite field: compute the lqnr in the base prime field, then use Tonelli-Shanks or Cipolla-Lehmer algorithm to compute square roots in the extension field.

Denote by SoS(q, k) this algorithm decomposing k as a sum of squares within any finite field F q . This is not always possible over infinite fields, but there Algorithm 19 still works anyway for the special case k = -1: just run it in the prime sub-field, since -1 must be in it.

Finite fields with even extension and φ

(A) = A H With φ (A) = A H , we need a matrix Y such that Y • Y H = Y • Y = -I .
This is not possible anymore over the complex field, but works for any even extension field, thanks to Algorithm 19. To see this, we consider next the finite field F q 2 , where q is a power of an arbitrary prime. Given a ∈ F q 2 , we adopt the convention that conjugation is given by the Frobenius automorphism:

a = a q . ( 17
)
The bar operator is F q -linear and has order 2 on

F q 2 . First, if -1 is a square in F q , then Y = √ -1 • I n works in F q 2 since then √ -1 = √ -1: Y • Y = √ -1 • I n √ -1 • I n = -I n .
Second, otherwise, q ≡ 3 mod 4 and then there exists a square root i of -1 in F q 2 , from Proposition 16. Further, one can build (a, b), both in the base field F q , such that

a 2 + b 2 = -1, from Algorithm 19. Finally Y = (a + ib) • I n in F q 2 n×n
is skew-unitary: indeed, since q ≡ 3 mod 4, we have that i q = i 3+4k = i 3 (-1) 2k = -i and, therefore, a

+ ib = (a + ib) q = a -ib. Finally Y • Y = (a + ib)(a -ib) • I n = -I n .

Any field with positive characteristic and φ (A) = A H

If -1 is a square in the base field, or within an even extension we have seen in Section 4.4 that there exists diagonal skew-unitary matrices. Otherwise, one can always resort to tridiagonal ones as in Section 4.3. For this, one can always build (a, b) in the base field such that a 2 + b 2 = -1 using Proposition 17. Then, Y = a b -b a ⊗ I n is a skew-unitary matrix. Indeed, since a and b live in the base field, they are invariant by the Frobenius automorphism. Therefore,

Y = a b -b a ⊗ I n = a -b b a ⊗ I n and Y • Y = (a 2 + b 2 ) • I n = -I n .
5 Towards a minimality result on the number of multiplications

Our Algorithm 12 computes the product of a matrix over a ring by its (hermitian) transpose using only 5 block multiplications and the (hermitian) transpose of one of these block multiplications. Here, we use consider some vector-spaces and thus, restrict ourselves to consider matrices over a field. We reformulate in this section the method introduced by de Groote in [START_REF] Hans Friedich Groote | On varieties of optimal algorithms for the computation of bilinear mappings II. Optimal algorithms for 2 × 2-matrix multiplication[END_REF] in order to prove that the tensor rank of the 2 × 2 matrix product is 7. This method is used to prove the following result: Theorem 23. There is no algorithm derived from non-commutative block 2 × 2 matrix product algorithms that computes the product of a matrix over a field by its (hermitian) transpose using only 4 block multiplications and the (hermitian) transpose of one of these block multiplications. This result does not state that it is never possible to multiply by the adjoint using fewer than 5 multiplications as shown by the following remark.

Remark 24. Over any ring with a square root i of -1, there is a computational scheme requiring 4 multiplications and computing the product of a 2 × 2-matrix by its transpose:

a b c d ( a c b d ) = (a+ib)(a-ib) × ac+bd (c+id)(c-id) = a 2 +b 2 × ac+bd c 2 +d 2 . ( 18
)
This is the case for instance over F 2 , where i = 1, or over the complex numbers. As this scheme requires for instance that aib = iba, at least some commutativity is required, thus in general it does not apply to block matrices and it is therefore not in the scope of Theorem 23.

The following section is devoted to shortly present the framework used in this part of our work.

The framework of bilinear maps encoded by tensors

We present de Groote's proof using a tensorial presentation of bilinear maps; we recall briefly this standpoint through the following well-known example of seven multiplications and we refer to [START_REF] Landsberg | Geometry and complexity theory[END_REF] for a complete introduction to this framework.

Example 25. Considered as 2 × 2 matrices, the matrix product C = A • B could be computed using Strassen algorithm by performing the following computations (see [START_REF] Strassen | Gaussian elimination is not optimal[END_REF]):

ρ 1 ← a 11 (b 12 -b 22 ), ρ 2 ← (a 11 + a 12 )b 22 , ρ 4 ← (a 12 -a 22 )(b 21 + b 22 ), ρ 3 ← (a 21 + a 22 )b 11 , ρ 5 ← (a 11 + a 22 )(b 11 + b 22 ), ρ 6 ← a 22 (b 21 -b 11 ), ρ 7 ← (a 21 -a 11 )(b 11 + b 12 ), c 11 c 12 c 21 c 22 = ρ 5 + ρ 4 -ρ 2 + ρ 6 ρ 6 + ρ 3 ρ 2 + ρ 1 ρ 5 + ρ 7 + ρ 1 -ρ 3 . (19)
With m, n, p equal to 2, this algorithm encodes a bilinear map:

F m×n × F n×p → F m×p , (A, B) → A • B. ( 20
)
We keep the indices m, n, p in this section for the sake of clarity in order to distinguish the different spaces involved in the sequel. The spaces F •×• can be endowed with the Frobenius product M, N = Trace(M • N ) that establishes an isomorphism between F •×• and its dual space F •×• ; hence, it allows for example to associate the trilinear form Trace(C • A • B) and the matrix multiplication [START_REF] Wedeniwski | Primality tests on commutator curves[END_REF]:

S | 3 : F m×n × F n×p × (F m×p ) → F, (A, B, C ) → C, A • B . (21)
As by construction, the space of trilinear forms is the canonical dual space of order three tensor products, we could encode the Strassen multiplication algorithm [START_REF] Strassen | Gaussian elimination is not optimal[END_REF] as the tensor S defined by:

S = 7 i=1 Σ i 1 ⊗Σ i 2 ⊗S 3 i = Σ i 1 ⊗Σ i 2 ⊗S 3 i = ( 1 0 0 0 )⊗ 0 1 0 -1 ⊗( 0 0 1 1 )+ ( 1 1 0 0 )⊗( 0 0 0 1 )⊗ -1 0 1 0 +( 0 0 1 1 )⊗( 1 0 0 0 )⊗ 0 1 0 -1 + 0 1 0 -1 ⊗( 0 0 1 1 )⊗( 1 0 0 0 )+ ( 1 0 0 1 )⊗( 1 0 0 1 )⊗( 1 0 0 1 )+( 0 0 0 1 )⊗ -1 0 1 0 ⊗( 1 1 0 0 )+ -1 0 1 0 ⊗( 1 1 0 0 )⊗( 0 0 0 1 ) (22)
in

(F m×n ) ⊗ (F n×p ) ⊗ F m×p with m = n = p = 2.
Remark that-as introduced in the above Equation ( 22)-we are going to use in the sequel the Einstein summation convention in order to simplify the forthcoming notations (according to this convention, when an index variable appears twice in a term and is not otherwise defined, it represents in fact the sum of that term over all the values of the index).

Starting from the tensor representation S of our algorithm, we could consider several contractions that are the main objects manipulated in the sequel.

Flattening tensors and isotropies

The complete contraction S | 3 (A ⊗ B ⊗ C ) is defined as the following map:

((F m×n ) ⊗ (F n×p ) ⊗ F m×p ) ⊗ (F m×n ⊗ F n×p ⊗ (F m×p ) ) → F, Σ i 1 ⊗Σ i 2 ⊗S 3 i ⊗ (A ⊗ B ⊗ C ) → Σ i 1 , A Σ i 2 , B S 3 i , C . (23)
We already saw informally in the previous section that this complete contraction is Trace(A • B • C) and we recall in the following remark some of its basic properties.

Remark 26. Given three invertible matrices:

α ∈ F m×m , β ∈ F p×p , γ ∈ F n×n (24)
that encodes changes of basis, the trace Trace(A • B • C) is equal to:

Trace (A • B • C) = Trace(C • A • B) = Trace(B • C • A), and Trace α -1 • A • β • β -1 • B • γ • γ -1 • C • α . ( 25
)
These relations illustrate the following theorem:

Theorem 27 ([12, § 2.8]). The isotropy group of the n × n matrix multiplication tensor is psl

± (F n ) ×3 S 3
, where psl stands for the group of matrices of determinant ±1 and S 3 for the symmetric group on 3 elements.

The following classical statement redefines the sandwiching isotropy on a matrix multiplication tensor:

Definition 28. Given g = (α × β × γ) in psl ± (F n ) ×3
, its action g S on a tensor S is given by g

(Σ i 1 ⊗ Σ i 2 ⊗ S 3 i )
where each summands is equal to:

(α -1 ) • Σ i 1 • β ⊗ (β -1 ) • Σ i 2 • γ ⊗ (γ -1 ) • S 3 i • α , ∀i fixed. (26)
These isotropies will be used later; for the moment, let us now focus our attention on the very specific standpoint on which is based the forthcoming developments: flattenings.

Definition 29. Given a tensor S, the third flattening S | 1 3 (a.k.a. third 1contraction) of the tensor S is:

S | 1 3 : F m×p → (F m×n ) ⊗ (F n×p ) , M → M, S 3 i Σ i 1 ⊗ Σ i 2 .
(

) 27
Example 30. To illustrate this definition and some important technicalities, let us consider Im S | 1 3 the image of the Strassen tensor [START_REF] Xianyi | an Optimized BLAS library[END_REF] flattening: this is a subspace of (F m×n ) ⊗ (F n×p ) . More precisely, let us first consider only the fifth summand in Equation [START_REF] Xianyi | an Optimized BLAS library[END_REF] and the image of its third flattening:

Im ( 1 0 0 1 )⊗( 1 0 0 1 )⊗( 1 0 0 1 ) | 1 3 = c11+c22 0 0 c11+c22 0 0 0 0 0 0 0 0 c11+c22 0 0 c11+c22 ∀(c 11 , c 22 ) ∈ F 2 . ( 28
)
The indeterminates c 11 and c 22 keep track of the domain of the flattening. The 4 × 4 right-hand side matrix in above expression should not be confused with the Kronecker product I 2×2 ⊗ I 2×2 involved in the left-hand side. In fact, the result I 4×4 of this Kronecker product is of classical matrix rank 4 while the rank in (F m×n ) ⊗ (F n×p ) of the right-hand side matrix (28) is 1 by construction. Hence, even if we present in this section the elements of F •×• as matrices, we use a vectorization (e.g. F m×n F mn ) of these matrices in order to perform correctly our computations. Taking this standpoint into account we obtain the following description of the whole Strassen third flattening image as:

Im S | 1 3 =     c 11 c 12 0 0 0 0 c 11 c 12 c 21 c 22 0 0 0 0 c 21 c 22     (29)
that could be guessed almost without computation. In fact, this right-hand side matrix is just the matrix of the bilinear form defining the trilinear encoding of the matrix product:

Trace(A • B • C ) = a 11 a 12 a 21 a 22     c 11 c 12 0 0 0 0 c 11 c 12 c 21 c 22 0 0 0 0 c 21 c 22         b 11 b 12 b 21 b 22     . (30)
Hence, the flattening Im S | 1 3 is a canonical description of the matrix product independent from the algorithm/tensor used to encode it; in particular, it is an invariant under the action of the isotropies introduced in Definition 28. We are going to use these objects and their properties in the following section.

Presentation of de Groote's method

We are interested in a situation where, given a bilinear map, a classical representation by a tensor C is known and we wish to disprove the existence of a tensor representation of a given rank. Inspired by Steinitz exchange theorem, de Groote introduced in [10, § 1.5] the following definition to handle this issue.

Definition 31. Given a tensor T encoding a bilinear map whose codomain is M and given q rank-one elements P i , linearly independent in (F m×n ) ⊗ (F n×p ) , let us denote by L(Im T | 1 3 , P 1 , . . . , P q ) the linear subspace of M defined by:

LinearSpan M ∈ M | ∃(u 1 , . . . , u q ) ∈ F q , Rank (F m×n ) ⊗(F n×p ) T | 1 3 (M ) + u j P j = 1 . ( 31
)
We introduced the notation M for the codomain of the considered bilinear map in order to highlight the fact that it is isomorphic-via the Frobenius isomorphism-to the domain of the flattening and to show how it is used in the following.

The following proposition allows to construct an effective test that checks if there exists a tensor of rank dim M + q that defines the considered bilinear map.

Proposition 32. If there exists a tensor T of rank dim M + q encoding a bilinear map with codomain M then there are q rank-one elements P i linearly independent in (F m×n ) ⊗ (F n×p ) such that L(Im T | 1 3 , P 1 , . . . , P q ) is M. Proof. Let us assume that there exists a tensor C encoding the considered bilinear map, that its tensor rank ρ is greater than dim M + q and that it is defined by the sum Q i ⊗ R i . Remark that the set (R i ) 1≤i≤ρ is a generating set of the space M by hypothesis.

Suppose now that there exists a tensor T of rank r = dim M + q encoding with fewer summands the same considered bilinear map:

T = P k ⊗ S k , S k ∈ M ⊂ F m×p , P k ∈ (F m×n ) ⊗ (F n×p ) . ( 32
)
The elements P k are linearly independent, otherwise T could be expressed with even fewer terms. Furthermore, there is a subset of (S i ) 1≤i≤dim M that is a base of M (otherwise, T could not encode the same bilinear map as C).

Suppose that we reorder this set so that the base is (S i+q ) 1≤i≤dim M . By invariance of the flattening map, the following relations hold:

∀ M ∈ M, C | 1 3 (M ) = M, R k Q k = T | 1 3 (M ) = M, S k P k . (33)
By introducing a base (B i ) 1≤i≤dim M of M, we could summarize this situation under a matricial standpoint as follows:

          P 1 . . . P q B 1 , R k Q k . . . B dim M , R k Q k           = I q×q 0 C D           P 1 . . . P q P q+1 . . . P q+dim M           (34)
where the matrices C and D are such that:

∀ i ∈ {1, . . . , dim M}, C ij = B i , S j , ∀ j ∈ {1, . . . , q}, D ij = B i , S q+j , ∀ j ∈ {1, . . . , dim M}. (35)
As, by hypothesis, (S i+q ) 1≤i≤dim M is a basis of M, the matrix D is invertible and we could rewrite Equation (34) as follows:

U = -D -1 • C, V = D -1 , I q×q 0 U V           P 1 . . . P q B 1 , R k Q k . . . B dim M , R k Q k           =           P 1 . . . P q P q+1 . . . P q+dim M           . (36)
The dim M lines of the (U, V ) matrices in Equation (36) give us dim M vectors (u j 1 , . . . , u j q , v j 1 , . . . , v j dim M ) such that for all j in {1, . . . , dim M}

Rank (F m×n ) ⊗(F n×p ) u j i P i + v j i B i , R k Q k = Rank (F m×n ) ⊗(F n×p ) P q+j = 1. ( 37
)
To conclude, we remark that these last relations show that all the matrices v j i B i are in the subspace L(Im T | 1 3 , P 1 , . . . , P q ). As they are dim M independent linear combinations of basis elements of M, these matrices form another of its bases and thus the subspace L(Im T | 1 3 , P 1 , . . . , P q ) is equal to M.

Adaptation to the Hermitian case

In order to use Proposition 32 to prove Theorem 23, we have to show that for any element P = P ⊗ Q in (F m×n ) ⊗ (F n×p ) the subspace L(Im T | 1 3 , P, P H ) is not equal to M (with P H = Q H ⊗ P H ). This vector-space M is a 3 dimensional vector-space spanned by all the outputs of our bilinear map. Let us start to make this strategy more precise by the following remark.

Remark 33. A classical block version of bilinear algorithm (e.g. [6, § 6.3.1]) computing the product of a matrix by its adjoint is:

A11 A12 A21 A22 A11 H A21 H A12 H A22 H = A11A11 H +A12A12 H × A21A11 H +A22A12 H A21A21 H +A22A22 H . ( 38
)
As the result of this algorithm is self-adjoint, by Lemma 9 there is no need to compute the top-right coefficient and thus, we conclude that the dimension of M is at most 3. Hence, there exists a bilinear map encoded by a tensor H = Σ i 1 ⊗ Σ i 2 ⊗ S 3 i of rank 6 that computes the product of a matrix by its hermitian transpose and the image of its third flattening H | 1 3 is

Im H | 1 3 (F m×p ) = c1 0 0 0 0 0 c1 0 c2 c3 0 0 0 0 c2 c3 . (39)
We need a last standard definition in order to classify all possible tensors P considered in the sequel.

Definition 34. Given a tensor P decomposable as sum of rank-one tensors:

P = q i=1 ⊗ s j=1 P ij where P ij are matrices. ( 40
)
The list [(Rank P ij ) j=1...s ] i=1...q is called the type of tensor P.

Remark 35. In our situation q is one, s is two and the P ij are 2 × 2 matrices; hence, the tensor P could only have type

[(1, 1)], [(1, 2)], [(2, 1)] or [(2, 2)].
We also use the isotropies presented in Definition 28 in order to simplify as much as possible the tensor P as illustrated in the proof of the following statement. Furthermore, let us first introduce several notations:

• we denote by G x=0 the matrix G in which the indeterminate x is replaced by 0;

• the determinant of the matrix

Gi,j G i,l G k,j G k,l is denoted by Det [i,j|k,l] (G);
• as the rank of a matrix is invariant under elementary row and columns operations, we also use the notation G|[i, j| ] for the matrix resulting from the addition to the ith line of G of its jth line multiplied by .

Lemma 36. There is no tensor P in (F m×n ) ⊗ (F n×p ) of type (1, i) with i equals to 1 or 2 such that the subspace L(H | 1 3 , P, P H ) is equal to M.

Proof. Let us consider a tensor P = A ⊗ B of type (1, i) with i = 1, 2. As the first component A is of rank one, there exists two vectors such that:

A = a 1 a 2 ⊗ b 1 b 2 = a 1 b 1 a 2 b 1 a 1 b 2 a 2 b 2 . ( 41
)
If the coefficient a 2 is zero, we choose a matrix β as the identity matrix. If the coefficient a 1 is zero, we could consider a permutation matrix β = ( 0 1 1 0 ); otherwise, if s = a 1 a 1 + a 2 a 2 = 0, consider:

β = a 1 a 2 a 2 -a 1 . ( 42
)
Then we have both

β • β H = s • I 2×2 and A • β = sb 1 0 sb 2 0 .
Hence, in any of these cases, there always exists a matrix β, which inverse is a multiple of its hermitian transpose, such that the isotropy g defined by (I m×m × β × I n×n ) satisfies the following properties:

g P = u 1 0 u 2 0 ⊗ v 1 v 2 v 3 v 4 , g P H = 1 s (g P) H (43)
With the above notations, conventions and isotropy's action, given any M in M, the 4 × 4 matrix H | 1 3 (M ) + y 1 P + y 2 P H is:

    y1u1v1 + y2v1 H u1 H + c1 y1u1v2 + y2v1 H u2 H y1u1v3 y1u1v4 y2v3 H u1 H y2v3 H u2 H c1 0 y1u2v1 + y2v2 H u1 H + c2 y1u2v2 + y2v2 H u2 H + c3 y1u2v3 y1u2v4 y2v4 H u1 H y2v4 H u2 H c2 c3     . ( 44
)
This matrix is supposed to be of rank one. Thus, all its 2 × 2 minors are equal to 0. Then, either c 1 or c 3 is equal to 0 and for any such P the subspace L(H | 1 3 , P, P H ) is thus not M.

There remains the case s = 0. Then let

β = a 1 -1 -a 2 0 a 1 . ( 45
)
The inverse of β is no longer related to β H anymore but β still transforms A into a single non-zero column matrix:

A • β = b 1 0 b 2 0
. Thus, for this β, the action of the isotropy g = (I m×m × β × I n×n ) is:

g P = b 1 0 b 2 0 ⊗ z 11 z 12 z 21 z 22 and g P H = U ⊗ V. ( 46
)
With the above notations, conventions and isotropy's action, given any M in M, the 4 × 4 matrix G = H | 1 3 (M ) + y 1 P + y 2 P H is thus:

y1b1z11+y2u11v11+c1 y1b1z12+y2u11v12 y1b1z21+y2u11v21 y1b1z22+y2u11v22 y2u12v11 y2u12v12 y2u12v21+c1 y2u12v22 y1b2z11+y2u21v11+c2 y1b2z12+y2u21v12+c3 y1b2z21+y2u21v21 y1b2z22+y2u21v22 y2u22v11 y2u22v12 y2u22v21+c2 y2u22v22+c3 . ( 47
)
This matrix is supposed to be of rank one in L(H | 1 3 , P, P H ). Thus, all its 2 × 2 minors are equal to 0.

On the one hand, if y 2 is zero, then the constraints of Equation ( 47) show that Det [2,3|4,4] (G y2=0 ) = c 1 c 3 is equal to 0. On the other hand, if y 2 is different from 0 then the minor Det [2,2|4,4] (G) is equal to c 3 y 2 u 12 v 12 and supposed to be equal to zero by hypothesis. We also have that the minor Det [2,1|4,4] (G) is equal to c 3 y 2 u 12 v 11 and is also supposed to be equal to zero by hypothesis. We are going to explore all the consequences induced by this constraint.

u12 = 0, v12 = 0 → c3 = 0, u12 = 0, v12 = 0 → Det [2,1|4,3] (Gu 12 =0,v 12 =0) = y2u22v11c1 = 0, u12 = 0, v12 = 0, u22 = 0 → Det [2,3|4,4] (Gu 12 =0,v 12 =0,u 22 =0) = c1c3 = 0, u12 = 0, v12 = 0, v11 = 0 see thereafter u12 = 0, v12 = 0 → Det [2,2|4,3] (Gu 12 =0) = -y2v12u22c1 = 0, u12 = 0, v12 = 0, u22 = 0 → c1 = 0, u12 = 0, v12 = 0, u22 = 0 → Det [2,3|4,4] (Gu 12 =0,u 22 =0) = c1c3 = 0. (48)
Now, if u 12 is different from 0, then from the first two minors, either The relations v 12 = v 11 = 0 hold or c 3 = 0.

Further, not both b 1 and b 2 can be zero, otherwise the tensor is of rank 0. W.l.o.g., suppose that b 2 = 0 and let

G = G v12=0,v11=0 |[1, 3|-b 1 /b 2 ]. Then Det [1,2|2,4] (G ) = (-b 1 /b 2 )y 2 c 3 u 12 v 22
and either b 1 = 0 or v 22 = 0 or c 3 = 0. This gives the following distinctions (recall that now u 12 = 0 and y 2 = 0):

b 1 = 0, → Det [1,1|2,4] (G b1=0 ) = y 2 u 12 v 22 c 1 = 0, b 1 = 0, v 22 = 0 → c 1 = 0, b 1 = 0, v 22 = 0 → Det [1,1|4,4] (G ) = c 1 c 3 = 0, b 1 = 0, v 22 = 0 → Det [1,2|4,4] (G ) = (-b 1 /b 2 )c 3 2 = 0, (49)
approach for product of matrix by its adjoint yields a 2M method using only 2 real products:

Algorithm 38 2M multiplication in an extension

Input: A commutative ring K and one of its extensions E; Input: A ∈ K m×n and B ∈ K m×n ; Input: φ an involutive matrix antihomomorphism of E.

Input: i ∈ E, commuting with K, and such that = iφ(i) ∈ K; Output:

(A + iB) • φ (A + iB) ∈ E m×m . Let H = A • φ (B) ∈ K m×m ; Let G = (A + B) • φ (A + B) ∈ K m×m ; return (G -H -φ (H)) + Hφ(i) + iφ (H).
Lemma 39. Algorithm 38 is correct. It costs 2MM K ω (n)+o(n ω ) operations over the base ring K. 2b) and (2c). As i commutes with K, we also have that

Proof. Let M = (A + iB) • φ (A + iB). By Lemma 4, we have that φ(iB) = φ (B)φ(i) = φ(i)φ (B). Thus, M = A • φ (A) + A • φ (B)φ(i) + iB • φ (A) + iB • φ (B)φ(i), by Equations (
M = A • φ (A) + B • φ (B) + A • φ (B)φ(i) + iB • φ (A). By Equations (2a) and (2c), we have that φ (H) = φ (φ (B)) • φ (A) = B • φ (A). Finally, G = A • φ (A) + A • φ (B) + B • φ (A) + B • φ (B) as φ( ) = φ(φ(i))φ(i) = iφ(i) = . Therefore G -H -φ (H) = G -H -φ (H) = A • φ (A) + B • φ (B) and M = G + Hφ(i) + iφ (H).
Example 40. For instance, if K = R, E = C, then i = √ -1 satisfies the conditions of Algorithm 38 for both cases when φ is the transposition or the conjugate transposition. Therefore, we obtain the multiplications of a matrix by its adjoint, whether it be the transpose or the conjugate transpose, in 2MM R ω + o(n ω ) operations in R. The classical divide and conquer algorithm, see e.g. [6, § 6.3.1], works directly over C and uses the equivalent of 2 2 ω -4 complex floating point n × n matrix products. Using the 3M method for the complex products, this algorithm uses overall 6 2 ω -4 MM R ω (n) + o(n ω ) operations in R. Finally, Algorithm 12 costs 2 2 ω -3 complex multiplications for a leading term bounded by 6 2 ω -3 MM R ω (n), improving over 2MM R ω for ω > log 2 (6) ≈ 2.585, but this does not apply to the conjugate transpose case. This is summarized in Table 2, also replacing ω by 3 or log 2 [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF] to illustrate the situation for the main feasible exponents.

The quaternion algebra

Given a field K of characteristic not 2, the K-algebra of quaternions H(K) is the K-vector space of all formal linear combinations: only when ω < log 2 (7) ≈ 2.81 and 2 < 6 2 ω -3 only when ω < log 2 (6) ≈ 2.585.

x 1 + x 2 i + x 3 j + x 4 k, (x 1 , . . . , x 4 ) ∈ K 4 , (55) Problem Alg
. MM 3 (n) MM log 2 7 (n) MM ω (n) A • B naive 8n 3 4 MM R log 2 (7) (n) 4 MM R ω (n) 3M 6n 3 3 MM R log 2 (7) (n) 3 MM R ω (n) A • A H Alg. 38 4n 3 2 MM R log 2 (7) (n) 2 MM R ω (n) [6] 3n 3 2 MM R log 2 (7) (n) 6 2 ω -4 MM R ω (n) A • A Alg. 38 4n 3 2 MM R log 2 (7) (n) 2 MM R ω (n) [6] 3n 3 2 MM R log 2 (7) (n) 6 2 ω -4 MM R ω (n) Alg. 12 2.4n 3 3 2 MM R log 2 (7) (n) 6 2 ω -3 MM R ω (n)
the non-commutative multiplication being defined by the bilinear extensions of the relations:

i 2 = j 2 = k 2 = ijk = -1. (56)
The quaternions can also be seen as a degree 2 extension of a degree 2 extension, but a non-commutative one. Therefore the 3M or §2M techniques of Section 6.1 only apply directly for the first degree 2 extension, while the second extension would require 4 multiplications. This gives 4 × 3 = 12 (resp. 4 × 2 = 8) multiplications in the base field for a general matrix multiplication (resp. a multiplication of a matrix by its transpose or conjugate transpose). For the former case, there exist actually algorithms using only 8 multiplications instead of 12. For the latter case, we present algorithms using only 7 multiplications for the transpose case and only 6 multiplications for the conjugate transpose case instead of 8.

Quaternions' multiplication

The multiplication of quaternions is (x 1 +x 2 i+x 3 j+x 4 k)(y 1 +y 2 i+y 3 j+y 4 k) = (w 1 + w 2 i + w 3 j + w 4 k), with:

w 1 = x 1 y 1 -x 2 y 2 -x 3 y 3 -x 4 y 4 (57)
w 2 = x 1 y 2 + x 2 y 1 + x 3 y 4 -x 4 y 3 (58)
w 3 = x 1 y 3 -x 2 y 4 + x 3 y 1 + x 4 y 2 (59)
w 4 = x 1 y 4 + x 2 y 3 -x 3 y 2 + x 4 y 1 (60)
Fiduccia showed in [START_REF] Fiduccia | Fast matrix multiplication[END_REF] how to compute this product with only 10 field multiplications and 25 additions, cleverly using Gauß' trick for the multiplication of complex numbers in three multiplications. Regarding the minimal number of base field operation required for the multiplication of quaternions, de Groote shows in [START_REF] Friedrich Groote | On the complexity of quaternion multiplication[END_REF] that 10 multiplications is minimal to compute both X • Y and Y • X. In addition, over the reals and the rationals, the minimal number of multiplications is 8 [START_REF] Howell | The complexity of the quaternion product[END_REF] and [START_REF] Hans Friedich Groote | On varieties of optimal algorithms for the computation of bilinear mappings II. Optimal algorithms for 2 × 2-matrix multiplication[END_REF]Proposition 1.7]. The algorithm of [START_REF] Howell | The complexity of the quaternion product[END_REF], requiring also 28 additions, is recalled in Algorithm 41.

Algorithm 41 Howell-Lafon quaternion multiplication

Input: x = x 1 + x 2 i + x 3 j + x 4 k ∈ H(K), y = y 1 + y 2 i + y 3 j + y 4 k ∈ H(K) Output: xy ∈ H(K) Q 1 = (x 1 + x 2 )(y 1 + y 2 ), Q 2 = (x 4 -x 3 )(y 3 -y 4 ) Q 3 = (x 2 -x 1 )(y 3 + y 4 ), Q 4 = (x 3 + x 4 )(y 2 -y 1 ) Q 5 = (x 2 + x 4 )(y 2 + y 3 ), Q 6 = (x 2 -x 4 )(y 2 -y 3 ) Q 7 = (x 1 + x 3 )(y 1 -y 4 ), Q 8 = (x 1 -x 3 )(y 1 + y 4 ) T 1 = Q 5 + Q 6 , T 2 = Q 7 + Q 8 T 3 = Q 5 -Q 6 , T 4 = Q 7 -Q 8 T 5 = T 2 -T 1 , T 6 = T 1 + T 2 T 7 = T 3 + T 4 , T 8 = T 3 -T 4 w 1 = Q 2 + T 5 /2, w 2 = Q 1 -T 6 /2 w 3 = T 7 /2 -Q 3 , w 4 = T 8 /2 -Q 4 return w 1 + w 2 i + w 3 j + w 4 k.
Proposition 42. Algorithm 41 extends to the case of quaternions with matrix coefficients. If matrix multiplication over the base field costs MM K ω (n) field operations for n×n matrices and the field matrix addition O(n 2 ) field additions, then the dominant cost of Algorithm 41 applied to matrices is bounded by

8MM K ω (n).
Proof. Correctness is by inspection since 2 is invertible in K of characteristic different from 2. The complexity bound is just the fact that the Algorithm performs 8 multiplications of matrices with coefficients in the base field.

The lowest number of multiplications required to multiply two quaternions being 8, Proposition 42 is the best possible result while keeping the view of the matrices as two quaternions with base field matrix coefficients, X 1 , X 2 , X 3 , X 4 and Y 1 , Y 2 , Y 3 , Y 4 . The alternative is to use a matrix with quaternion coefficients and use classical fast matrix algorithms. Next, we see the different alternatives for the multiplication by an adjoint.

Multiplication of a quaternion matrix by its transpose

We now propose several methods to multiply a quaternion matrix by its transpose:

1. First a "7M" method which considers a quaternion with matrix coefficients, and thus reduces everything to seven general matrix multiplications over the base field.

2. Second, one can consider a matrix with quaternion coefficients and just apply any matrix multiplication algorithm where multiplication of coefficients is that of Algorithm 41.

7M method: a quaternion with matrix coefficients. Many simplifications used in computing the square of a quaternion no longer apply when computing the product of a quaternion matrix by its transpose, due to noncommutativity of the matrix product. For A, B, C, D ∈ K m×n ,

(A + Bi + Cj + Dk)(A + B i + C j + D k) = S 1 + S 2 i + S 3 j + S 4 k (61)
where:

S 1 = AA -BB -CC -DD (62)
S 2 = (AB + BA ) + (CD -DC ) (63)
S 3 = (AC + CA ) + (DB -BD ) (64)
S 4 = (AD + DA ) + (BC -CB ) (65)
Using Munro's trick twice, this can be computed with 7 multiplications over the field K and 17 additions (6 of which are half-additions), as shown in Algorithm 43.

Open question 1. Multiply a quaternion with matrix coefficients by its transpose in fewer than 7 multiplications.

Using matrices of quaternions and divide and conquer. In the following, for the sake of simplicity, we will consider only square matrices. Here, we consider instead a matrix with quaternion coefficients and perform a matrixmatrix product: since the quaternions are not commutative, then M M is not necessarily symmetric, one has to compute both the top right and bottom left corners of the product. Thus no gain is obvious between computing M M and M N that way and Algorithm 12 is a priori useless in this case, as remarked in the (counter)-Examples 6. The idea is thus to use a non symmetric algorithm, applied to M and M . The baseline cost would then again be MM H(K) symmetric products instead of 2 for the algorithms of [START_REF] Strassen | Gaussian elimination is not optimal[END_REF][START_REF] Winograd | La complexité des calculs numériques[END_REF]. In the quaternion case, as transposition is not an antihomomorphism, one cannot use the symmetries directly to save computations as in general (M • N ) = N • M . But if one is willing to recompute N • M then the algorithm still works. This yields to Algorithm 44 which requires 7 multiplications instead of 5, 15 additions instead of 9, and 4 multiplications by Y or Y .

ω (n) = 8MM K ω (n).
Algorithm 44 Product of a matrix by its non antihomomorphic transpose Input: A ∈ R m×n (with even m and n for the sake of simplicity);

Input: Y ∈ R n 2 × n 2 such that Y • Y = -I n 2 ; Output: A • A . Split A = A11 A12 A21 A22 where A 11 is in R m 2 × n 2 S 1 ← (A 21 -A 11 ) • Y ST 1 ← Y • (A 21 -A 11 ) S 2 ← A 22 -A 21 • Y ST 2 ← A 22 -Y • A 21 S 3 ← S 1 -A 22 ST 3 ← ST 1 -A 22 S 4 ← S 3 + A 12 ST 4 ← ST 3 + A 12 P 1 ← A 11 • A 11 P 2 ← A 12 • A 12 P 3 ← A 22 • ST 4 P T 3 ← S 4 • A 22 P 4 ← S 1 • ST 2 P T 4 ← S 2 • ST 1 P 5 ← S 3 • ST 3 U 1 ← P 1 + P 5 U 3 ← P 1 + P 2 U 2 ← U 1 + P 4 U T 2 ← U 1 + P T 4 U 4 ← U 2 + P 3 U T 4 ← U T 2 + P T 3 U 5 ← U 2 + P T 4 return U3 U T 4 U4 U5 .
Now, it turns out that transposition is still antihomomorphic if one of the matrices has its coefficients in the base field, as shown by Lemma 45.

Lemma 45. Let A be in H(K) m×k and Y in K k×n , then (A • Y ) = Y • A .
Proof. Since the coefficients of Y are in the base field, they commute with the quaternions. Therefore, ∀i, j, k a jk y ki = k y ki a jk . Now, Section 4.3 shows that for any quaternion algebra in positive characteristic, there exist a matrix Y , in the base field, such that Y • Y = -I n 2 . Therefore, in this case, Lemma 45 shows that in Algorithm 44, ST 1 = S 1 , ST 2 = S 2 , ST 3 = S 3 , ST 4 = S 4 . This shows that not only P 1 and P 2 are multiplications of a matrix by its transpose, but also P 5 = S 3 • S 3 . Finally, Algorithm 44 thus requires three recursive calls and four general multiplications. This is:

P (n) ≤ 3P ( n 2 ) + 4MM H(K) ω ( n 2 ) + o(n ω ) ≤ 3P ( n 2 ) + 32MM K ω ( n 2 ) + o(n ω ) (69)
and Equation ( 68) is modified as:

P (n) ≤ 32 2 ω -3 MM K ω (n) + o(n ω ) (70)
As expected this is again 8MM K log 2 (7) (n) if a Strassen-like algorithm is also used for the baseline over the field and ω = log 2 [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF]. This is again worse if ω < log 2 [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF], but better, and only (6

+ 2 5 )MM K 3 (n), if ω = 3.

Multiplication of a quaternion matrix by its adjoint

We now deal with the case of the product of a quaternion matrix with its conjugate transpose. This operator is now an antihomomorphism which allows us to save some computations as in Algorithm 12. Here also we distinguish the matrix of quaternions from the quaternion with matrix coefficients.

Scalar case. The quaternion conjugation satisfies XY = Y X. Therefore, we have:

(a + bi + cj + dk)(a + bi + cj + dk) = a 2 + b 2 + c 2 + d 2 (71)
For matrices again simplifications do not occur and the product is then more complex.

Using matrices of quaternions. Now M • M H is a hermitian matrix and Algorithm 12 works over H(K). For this, one needs to find a skew-unitary matrix in H(K). This is impossible in H(C), but always possible in the quaternions over fields of positive characteristic using sums of squares and Equation (71). Suppose we use a generic matrix multiplication algorithm over the quaternions with cost bound equivalent to MM H(K) ω (n) field operations. Then our Algorithm 12 can multiply a matrix of a quaternions by its conjugate transpose with a dominant complexity term bounded by Now for the quaternions, the best algorithm to multiply any two matrices of quaternions is given by Proposition 42 and uses 8MM K ω (n) field operations if the base field matrix multiplication uses MM K ω (n). We thus have proven:

Corollary 46. Algorithm 12 multiplies a quaternion matrix by its conjugate transpose with dominant cost bounded by Directly using quaternions with matrix coefficients. Using a quaternion with matrix coefficients over the field, we have:

M • M H = (A + Bi + Cj + Dk) • (A + Bi + Cj + Dk) H = (A + Bi + Cj + Dk)(A -B i -C j -D k) = H 1 + H 2 i + H 3 j + H 4 k (72)
where:

H 1 = AA + BB + CC + DD (73)
H 2 = (BA -AB ) + (DC -CD ) (74)
H 3 = (CA -AC ) + (BD -DB ) (75)
H 4 = (DA -AD ) + (CB -BC ) (76)
Note that H 1 is symmetric and H 2 , H 3 , H 4 are skew-symmetric.

The properties of the transpose in the field shows that these can be computed with 4 + 3 * 2 = 10 multiplications (including four squares).

Now consider E = A + Bi and F = C + Di, so that M = E + F j. This shows that:

M • M H = (E + F j) • (E H + jF H ) = (EE H + F F H ) + (F jE H -EjF H ) (77)
Then, we have:

F jE H = (C + Di)j(A -B i) = (CA -DB )j + (DA + CB )k = Xj + Y k (78) -EjF H = (A + Bi)j(-C + D i) = (-AC + BD )j -(BC + AD )k = -X j -Y k (79)
Using Equations ( 78) and (79), we thus have Algorithm 47 which uses only 6 multiplications (one of which is a square) and a total of 14 additions, 7 of them being half-additions (On the one hand, 3 multiplications and 5 additions for Equations ( 78) and (79) overall, then 2 half-additions for Low (H 3 ) and Low (H 4 ); on the other hand, 2 multiplications, 1 square, 3 additions and 9 half-additions for Low (H 1 ) and Low (H 2 )).

Proposition 48. Algorithm 47 multiplies a quaternion matrix by its conjugate transpose with cost equivalent to 5 + 2 2 ω -3 MM K ω (n) base field operations.

Proof. Algorithm 47 uses 6 multiplications, one of which, Q 6 , is the product of a matrix in K by its transpose.

Open question 2. Multiply a quaternion with matrix coefficients by its Hermitian transpose in fewer than 6 multiplications (including 1 square), or with more squares.

Algorithm 47 Fast quaternion matrix multiplications by its adjoint

Input: A, B, C, D ∈ K m×n Output: M • M H ∈ H(K) m×m , for M = A + Bi + Cj + Dk. U 1 = A + B, U 2 = C + D Q 1 = CA , Q 2 = DB Q 3 = U 2 U 1 Q 4 = AB , Q 5 = CD Q 6 = (U 1 + U 2 )(U 1 + U 2 ) T 1 = Q 1 -Q 2 , T 2 = Q 4 + Q 5 T 3 = (Q 1 + Q 2 ) -Q 3 Low (H 1 ) = Low (Q 6 -(Q 3 + Q 3 ) -(T 2 + T 2 )) Low (H 2 ) = Low (T 2 -T 2 ) Low (H 3 ) = Low (T 1 -T 1 ) Low (H 4 ) = Low (T 3 -T 3 ) return H 1 + H 2 i + H 3 j + H 4 k.
Comparison. We summarize the results of this section about quaternion matrices in Table 3.

Algorithm into practice

This section reports on an implementation of Algorithm 12 over a prime field, as it is a core ingredient of any such computation in positive characteristic or over Z[i] or Q[i]. In order to reduce the memory footprint and increase the data locality of the computation, we first need to identify a memory placement and a scheduling of the tasks minimizing the temporary allocations. We thus propose in Table 4 and Figure 1 a memory placement and schedule for the operation C ← A • A using no more extra storage than the unused upper triangular part of the result C.

The more general operation C ← αA • A + βC, is referred to as SYRK (Symmetric Rank k update) in the blas api. Table 5 and Figure 2 propose a schedule requiring only one additional n/2 × n/2 temporary storage.

These algorithms have been implemented as the fsyrk routine in the open source fflas-ffpack library for dense linear algebra over a finite field [8, from commit 0a91d61e].

Figure 3 compares the computation speed in effective Gfops (a normalization, defined as n 3 /(10 9 × time)) of this implementation over Z/131071Z with that of the double precision blas routines dsyrk, the classical cubic-time routine over a finite field (calling dsyrk and performing modular reductions on the result), and the classical divide and conquer algorithm [6, § 6. (7) (n) 16 MM K ω (n) [START_REF] Howell | The complexity of the quaternion product[END_REF] 16n 3 8 MM K log 2 (7) (n) 8 MM K ω (n)

A • A H Alg. 47 10.8n 3 5.5 MM K log 2 (7) (n)

2 2 ω -3 + 5 MM K ω (n)
Alg. 12 6.4n 3 4 MM K log 2 (7) (n) (7) (n) 7 MM K ω (n) Eq. (68) (13 + 1

3 )n 3 8 MM K log 2 (7) (n) with gcc-9.2 on an Intel skylake i7-6700 running a Debian gnu/Linux system (v5.2.17).

The slight overhead of performing the modular reductions is quickly compensated by the speed-up of the sub-cubic algorithm (the threshold for a first recursive call is near n = 2000). The classical divide and conquer approach also speeds up the classical algorithm, but starting from a larger threshold, and hence at a slower pace. Lastly, the speed is merely identical modulo 131041, where square roots of -1 exist, thus showing the limited overhead of the preconditioning by the matrix Y .

Perspective

We made progresses in order to prove that five non-commutative products are necessary for the computation of the product of a 2 × 2 matrix by its adjoint, by applying de Groote's method to this context. However, we only prove that there is no algorithm, derived from a bilinear one, which uses 4 products and the adjoint of one of them. The case where more adjoints of already computed products could be used need to be ruled out in a similar manner. More generally, the possible existence of algorithms not originating from a bilinear algorithm is an even more challenging question.

Over the algebra of quaternions, the natural generalization of Howell and Lafon's algorithm to matrix coefficients yields the 7M and 6M algorithms for the transpose and conjugate transpose respectively. The recursive 5 products algorithm is only usable for the conjugate transpose case in positive characteristic and costs, as expected, half the cost of a general quaternion matrix product for ω = log 2 7. Yet for ω < 2.536 and for the case of transposition the 6M and 7M algorithms perform best. The minimality of the number 6 of multiplications to multiply a quaternion by its conjugate is an open question, as for the minimality of the number 7 of multiplications to multiply a quaternion by its transpose. For these questions, de Groote's method could provide an answer.

We proposed several algorithms for the product of a matrix by its adjoint, each of which improves by a constant factor the best known costs, depending on the algebraic nature of the field of coefficients and on the underlying matrix exponent to be chosen. When implemented in practice the comparison may become even more complex, as other parameters, such as memory access pattern or vectorization will come into play. Our first experiments show that these constant factor improvements do have a practical impact.

Definition 7 .

 7 The image φ (A) of a matrix A by an antihomomorphism is called the adjoint of A. Definition 8. Let A ∈ R m×n , we denote respectively by Low (A) and Up (A) the m × n lower and upper triangular parts of A, namely the matrices L and U verifying

  then half of the non-zero elements x in the base field of size p satisfy x p-1 4

  p

8 :

 8 return (a, ac mod p) Proposition 20. Algorithm 19 is correct and, under the erh, runs in expected time O log 3 (p) .

2 2 ω

 2 -3 MM H(K) ω (n) operations, by Theorem 13.

Figure 2 :

 2 Figure 2: dag of the tasks and their memory location for the computation of C ← αA • A + βC presented in Table5.

  Figure 2: dag of the tasks and their memory location for the computation of C ← αA • A + βC presented in Table5.

Figure 3 :

 3 Figure 3: Speed of an implementation of Algorithm 12

Table 2 :

 2 Multiplication of a matrix by its adjoint or transpose over C: leading term of the cost in number of arithmetic operations over R. Note that 2 < 6 2

ω -4

  Another approach is to use a divide and conquer strategy at the higher level:• M 11 • M 21 , M 12 • M 22 , M 21 • M 11 , M22 • M 12 by the baseline algorithm; • and M 11 • M 11 , M 12 • M 12 , M 21 • M 21 , M 22 • M 22 by recursive calls.

	cut M into	M 11 M 12 M 21 M 22	, and compute:

  3.1]. The fflas-ffpack library is linked with Openblas [22, v0.3.6] and compiled

		Alg.	MM3(n)	MM log 2 7 (n)	MMω(n)
	A • B	naive	32n 3	16 MM K log 2

Table 3 :

 3 Matrix Multiplication over H(K) n×n : leading term of the cost in number of operations over K. Note that2 2 ω -3 + 5 < 16 2 ω -3 only when ω < log 2 (29) -log 2 (5) ≈ 2.536. = (A 21 -A 11 ) • Y C 21 9 U 1 = P 1 + P 5 C 12 2 S 2 = A 22 -A 21 • Y C 12 Up(U 1 ) = Low(U 1 ) C 12 3 P 4 = S 2 • S 1 C 22 10 U 2 = U 1 + P 4 C 12 4 S 3 = S 1 -A 22 C 21 11 U 4 = U 2 + P 3 C 21 5 P 5 = S 3 • S 3 C 12 12 U 5 = U 2 + P 4 C 22 6 S 4 = S 3 + A 12 C 11 13 P 2 = A 12 • A 12 C 12 7 P 3 = A 22 • S 4 C 21 14 U 3 = P 1 + P 2 C 11 8 P 1 = A 11 • A 11 C 11

	# operation	loc. # operation	loc.
	1 S 1

Table 4 :

 4 Memory placement and schedule of tasks to compute the lower triangular part of C ← A • A when k ≤ n. The block C 12 of the output matrix is the only temporary used.Figure 1: dag of the tasks and their memory location for the computation ofC ← A • A presented in Table 4. = (A 21 -A 11 ) • Y tmp P 1 = αA 11 • A 11 tmp S 2 = A 22 -A 21 • Y C 12 U 1 = P 1 + P 5 C 12 Up(C 11 ) = Low(C 22 ) C 11 Up(U 1 ) = Low(U 1 ) C 12 P 4 = αS 2 • S 1 C 22 U 2 = U 1 + P 4 C 12 S 3 = S 1 -A 22 tmp U 4 = U 2 + P 3 C 21 P 5 = αS 3 • S 3 C 12 U 5 = U 2 + P 4 +βUp(C 11 ) C 22 S 4 = S 3 + A 12 tmp P 2 = αA 12 • A 12 + βC 11 C 11 P 3 = αA 22 • S 4 + βC 21 C 21 U 3 = P 1 + P 2 C 11

	32

Table 5 :

 5 Memory placement and schedule of tasks to compute the lower triangular part of C ← αA • A + βC when k ≤ n. The block C 12 of the output matrix as well as an n/2 × n/2 block tmp are used as temporary storage.

	C 11	C 22	C 12	tmp	C 21
	Up(C 11 )		S 2	S 1
		P 4		S 3
			P 5	S 4
				P 1	P 3
			U 1
			U 2
		U 5			U 4
	P 2
	U 3

ω -3 MM K ω (n) base field operations.

There remains the case u 12 = 0, v 12 = 0 and v 11 = 0 in Equation (48). Here also, w.l.o.g., suppose that b 2 = 0 and let

Thus, in any cases, for any such P, the set L(H | 1 3 , P, P H ) is not M. Note that a computational way to see this, is to perform a Gröbner basis computation, directly from Equation (47): for instance over C this gives that the relation c 1 2 c 3 2 = 0 must hold and that the set is not the full codomain.

According to Remark 35, the above computations deal with half of the cases to consider. We remark that, mutatis mutandis, similar computations exclude also the existence of an algorithm where P is of type (2, 1). We could consider now the last case.

Lemma 37. There is no tensor P in (F m×n ) ⊗ (F n×p ) of type (2, 2) such that the subspace L(H | 1 3 , P, P H ) is equal to M.

Proof. First, let us consider a tensor P of type (2, 2). Thus, there exists β such that the action of the isotropy g = (I m×m × β × I n×n ) is:

With the above notations, conventions and isotropy's action, given any M in M, the 4 × 4 matrix G = H | 1 3 (M ) + y 1 P + y 2 P H is:

This matrix is supposed to be of rank one in L(H | 1 3 , P, P H ). Thus, all its 2 × 2 minors are equal to 0.

A Gröbner basis computation over C shows in that case that the relations c 1 2 c 3 = c 1 c 2 c 3 = c 1 c 3 2 = 0 hold and this is sufficient to conclude. Nevertheless, we present a proof that does not require such computations and is valid for any field.

On the one hand, if y 2 = 0, then the constraints of Equation ( 52) show for instance that both Det [2,1|3,3]

On the other hand, if y 2 = 0 then consider the minor Det [2,3|3,4] (G), which is equal to c 1 y 2 u 21 v 22 and supposed to be equal to zero by hypothesis. We are going to explore all the consequences induced by this constraint. First, if y 1 = 0, then we have:

If y 1 and y 2 are both non-zero, then, the minor Det [1,1|2,2] (G|[1, 4|-1]) is equal to y 2 c 1 u 12 v 12 and supposed to be equal to zero by hypothesis. We are going to explore all the consequences induced by this constraint. Let G = G|[1, 4|-1]:

Thus, in any cases, for any such P, the set L(H | 1 3 , P, P H ) is not M.

The Proposition 32, together with the computations done in Lemma 37 and in Lemma 36 are sufficient to conclude the proof of Theorem 23. [START_REF] Dumas | Dense linear algebra over prime fields[END_REF] The case of field extensions via matrix polynomial arithmetic

The cost comparison in Table 1 is for matrices over an arbitrary ring with skew unitary matrices. When the ring is an extension, the input of the problem is a polynomial matrix over the base ring. Following the traditional equivalence between polynomial matrices and matrix polynomials, leads to alternative ways to multiply the matrix by its transpose, considering the product of two polynomials with matrix coefficients. More specifically, we will focus on degree two extensions, and compare the costs in terms of number of operations over the base ring.

The 2M method

Over the field C of complex numbers, the 3M method (Karatsuba) for general matrix multiplication reduces the number of multiplications of real matrices from 4 to 3 [START_REF] Higham | Stability of a method for multiplying complex matrices with three real matrix multiplications[END_REF]: if MM R ω (n) is the cost of multiplying n × n matrices over R, then the 3M method costs 3MM R ω (n) + o(n ω ) operations over R. Adapting this Algorithm 43 Fast quaternion matrix multiplications by its transpose

The cost of this divide and conquer strategy is then:

By Lemma 11, we have that C(n) ≤ 32 2 ω -4 MM K ω (n) + o(n ω ), and this is never better that 8MM K ω (n) (but equal when ω = 3 as expected). So this is thus useless too.

But the same strategy can be used with a Strassen-like algorithm instead. Now such algorithms, for instance those of [START_REF] Strassen | Gaussian elimination is not optimal[END_REF][START_REF] Winograd | La complexité des calculs numériques[END_REF], when applied to M and M , use two recursive calls and five normal multiplications. This is:

By Lemma 11, we obtain that this is

As expected this is again 8MM K log 2 (7) (n) if a Strassen-like algorithm is also used for the baseline over the field and ω = log 2 [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF]. This is worse if ω < log 2 [START_REF] Dumas | On fast multiplication of a matrix by its transpose[END_REF], but better, and only (6 + 2

3 )MM K 3 (n), if ω = 3.

Positive characteristic quaternions and transposition. In this case, Algorithm 12 is not usable. It nonetheless has an interesting feature: it has 3